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SECTION  I 


RESEARCH  OBJECTIVES 


1.1  Overall  Objectives 

The  subject  of  combustion  instability  in  solid  propellant  rocket  motors 
has  been  studied  by  many  investigators  for  the  past  thirty  years.  Today,  stable 
motors  are  designed  routinely  and  there  seems  to  be  little  apparent  concern  . 
over  combustion  instability.  Looking  into  the  history  of  the  developmental  pro¬ 
cess,  however,  one  finds  that  more  than  half  of  all  the  motors  developed  have 
been  found  unstable.  If  new  designs  for  more  powerful  motors  are  proposed,  a 
prediction  as  to  stability  would  be  difficult  unless  trial  and  error  procedures 
are  repeated  with  costly  experiments. 

The  current  practice  for  the  prediction  of  combustion  instability  is  based 
on  crude  approximations.  In  most  of  the  design  calculations,  one-dimensional 
analyses  are  predominantly  used.  Seldom  are  the  mean  flow  calculations  per¬ 
formed  using  the  most  modern  technology  -  computational  fluid  dynamics.  In  the 
cavity  of  the  solid  propellant  rocket  motor  are  the  extremely  complicated  fluid 
mechanics  problems  -  compressible  viscous  flow,  vortex  shedding,  turbulent  bound¬ 
ary  layers,  particle  damping  (two  phase  flow),  etc.  In  addition,  oscillatory 
motions  are  prevalent,  with  acoustic  and  hydrodynamic  wave  oscillations  coupled 
toegther.  The  boundary  conditions  for  this  flowfield  are  the  unsteady  respons¬ 
es  ^Of  'the  flame  zone  distribution  of  field  variables  (velocity,  pressure,  den¬ 
sity,  temperature,  and  fuel  fractions).  The  system  may  be  linearly  or  nonlin- 
early  unstable.  Steep-fronted  waves,  erosive  burning,  and  high  amplitude  re¬ 
sponses  as  related  to  the  velocity-coupling  are  the  recent  subjects  of  interest 
which  are  not  fully  understood.  It  Is  clear  that  the  motor  consists  of  both 
flame  zone  and  cavity,  and  that  these  two  regions  should  not  be  separated  in 
the  analysis.  Obviously,  this  is  beyond  the  current  state  of  the  art,  but  it 
Is  toward  this  goal  that  our  overall  objective  must  be  intended. 

Analytical  or  numerical  models  should  be  developed  such  that  future  high 
performance  motors  are  designed  using  the  computerized  procedure.  These  models 
must  be  based  on  adequate  multi-dimensional  governing  equations  and  numerical 
methods.  The  results  should  also  be  verified  by  experimental  measurements. 

The  final  product  should  then  be  facilitated  by  an  interactive  computer  graph¬ 
ics  system  displaying,  for  example,  computed  waterfall  data  from  which  suitable 
design  decisions  can  be  made. 


1.2  Specific  Objectives 

The  overall  objective  stated  above  represents  our  ultimate  goal.  The  re¬ 
search  performed  during  the  two  year  period  consists  of  (1)  mean  flow  calcula¬ 
tions  and  interactions  of  unsteady  acoustic  and  vortical  oscillations  in  an 
axisymmetric  cylindrical  cavity,  (2)  particle  damping  effects  on  combustion  in¬ 
stability,  and  (3)  unsteady  response  of  the  burning  surface  in  solid  propellant 
combustion. 

It  is  anticipated  that  these  investigations  will  contribute  to  the  founda¬ 
tions  on  which  further  advancements  can  be  made  In  the  future. 
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SECTION  2 


MEAN  FLOW  CALCULATIONS  AND  INTERACTIONS  OF  UNSTEADY 
ACOUSTIC  AND  VORTICAL  OSCILLATIONS  IN 
AXI SYMMETRIC  CYLINDRICAL  CAVITY 


2 . 1  Summary 

The  flowfield,  such  as  occurs  in  solid  propellant  rocket  motors,  offers  a 
fertile  ground  for  fundamental  research  in  fluid  mechanics  and  heat  transfer. 
Combustion  induces  not  only  the  mean  flowfield,  but  also  acoustic  pressure  os¬ 
cillations  and  possibly  vortex  fluctuations  together  with  turbulent  shear 
boundary  layers.  Furthermore,  shock  waves  are  commonplace  in  most  instances. 
Obviously,  a  most  rigorous  analysis  taking  into  account  all  of  these  phenomena 
would  be  difficult,  if  not  impossible.  However,  with  the  advent  of  the  elec¬ 
tronic  computer  and  the  modem  technology  of  numerical  methods,  it  has  become 
feasible  to  resolve  hitherto  unsolved  problems. 

Despite  difficulties  in  analytical  and  numerical  solutions  to  the  complex 
physical  phenomena  in  a  rocket  motor  chamber,  many  researchers  have  contributed 
to  the  advancement  of  analysis  and  design  of  successful  rocket  motors.  A  large 
body  of  literature  exists  relative  to  this  subject,  the  study  of  which  has  been 
pioneered  by  Crocco  [1],  Cantrell  and  Hart  [2],  Culick  [3,4],  and  others. 
Flandro  and  Jacobs  [5],  among  others,  have  noted  that  vortex  shedding  may  lead 
to  an  instability  in  solid  propellant  rocket  motors.  It  Is  quite  possible 
that  high  speed  mean  flows  also  affect  the  stability  [6]  significantly. 

The  basic  mathematical  formulations  of  combustion  instability  were  con¬ 
tributed  by  Culick  [3,4].  Recently,  it  has  been  observed  in  both  full-scale 
firings  and  cold  flow  simulations  that  interactions  of  acoustic  and  hydrodynam- 
ical  (vortical)  Instability  can  be  significant  [7-13],  Although  It  can  be  ar-r 
gued  that  the  hydrodynamic  instability  may  not  occur  In  high  Reynolds  numbers, 
the  turbulent  shear  layer  instabilities  have  been  found  to  be  affected  by  vari¬ 
ous  combinations  of  Strouhal  numbers  and  Reynolds  numbers.  The  acoustic  field 
may  Interact  with  vortex  motions  known  as  the  "feedback"  resulting  in  the  vor¬ 
tex  generated  sound  [14].  Some  studies  [15,16]  Indicate  that  the  vortices  may 
undergo  "clipping",  a  phenomenon  corresponding  to  the  vortex  disruption.  It 
is  also  possible  that  lateral  periodic  motion  of  the  vortex  street  known  as 
"jitters"  may  lead  to  partial  or  complete  escape  of  the  vortices  T 1 5 , 16 ] .  Whe¬ 
ther  these  conditions  prevail  in  large  rocket  motors  in  which  flow  separations 
at  interface  restrictors  or  inhibitors  are  likely  to  produce  vortex  motion 
must  be  clarified.  No  simple  models,  such  as  hyperbolic  tangent  velocity  pro¬ 
file  for  the  shear  layer  [11,17]  and  temporal  or  spatial  growth  theories  [18, 
19],  appear  to  be  adequate  for  the  interactions  of  acoustic  and  vortical  os¬ 
cillations  in  a  rocket  chamber. 

In  the  previous  papers  [20-22],  finite  element  applications  to  the  com¬ 
bustion  instability  analysis  were  discussed.  Although  a  rigorous  mathematical 
formulation  of  the  stability  integral  was  presented,  the  mean  flow  calcula¬ 
tions  did  not  include  turbulent  flows.  Since  the  turbulent  flowfield  is  in¬ 
volved  in  shear  boundary  layers  and  vortex  motions,  it  is  intended  that  this 
subject  be  considered  in  the  mean  flow  calculations  and,  subsequently,  in  the 
stability  integral.  Shock  waves  will  not  be  included  in  the  present  report . 
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The  numerical  results  for  certain  combinations  of  acoustic  and  vortical 
frequencies  indicate  that  stability  boundaries  for  acoustics-coupled  hydro- 
dynamic  oscillations  are  somewhat  similar  to  the  classical  hydrodynamic  sta¬ 
bility  boundaries,  but  they  occur  in  the  form  of  multiple  islands.  The  tur¬ 
bulent  flowfield  appears  to  contribute  toward  instability,  and  this  trend  in¬ 
creases  with  larger  transition  angles  of  the  rocket  motor  cross-section.  De¬ 
tails  of  mathematical  formulations,  numerical  calculations,  and  example  prob¬ 
lems  are  presented  in  Reference  [23]  or  Appendix  1. 
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SECTION  3 


PARTICLE  DAMPING  EFFECTS  ON  COMBUSTION  INSTABILITY 


3 . 1  Summary 

The  study  of  acoustic  energy  losses  due  to  aluminum  particles  in  the  so¬ 
lid  propellant  rocket  motor  combustion  chamber  has  been  carried  out  by  a  num¬ 
ber  of  investigators.  Epstein  and  Charhart  [1]  studied  the  absorption  of  sound 
in  suspensions  of  non-interacting  inert  spherical  particles  and  uniform  temp¬ 
eratures.  The  validity  of  this  investigation  was  subsequently  substantiated 
by  other  researchers  [2-4].  In  a  rocket  motor,  the  particles  are  neither  spher¬ 
ical  nor  inert  and  subjected  to  nonuniform  temperature  distributions  [5]. 

Despite  the  extensive  research  on  the  subject  of  acoustic  energy  dissipa¬ 
tion  due  to  particle  damping  in  the  rocket  motor  [5-10],  calculations  of  the 
stability  integral  arising  from  particle  damping  are  limited  to  simple  one-di¬ 
mensional  cases. 

Thus,  the  purpose  of  the  present  study  is  to  demonstrate  the  feasibility 
of  mathematical  formulations  and  numerical  calculations  via  finite  elements. 
Furthermore,  interactions  of  particle  damping  with  fluid  viscosity  and  heat 
transfer  are  included.  It  is  shown  that  additional  boundary  and  domain  terms 
arise  from  integrating  by  parts  "twice"  of  the  Green  function  stability  inte¬ 
gral  containing  the  momentum  equation. 

Simple  example  problems  of  two-dimensional  axisymmetric  geometries  are 
solved  and  compared  with  one-dimensional  approximations.  It  is  shown  that 
there  is  a  trend  toward  decrease  in  stability  with  a  decrease  in  frequencies. 
However,  the  one-dimensional  calculation  overestimates  the  stability  at  high¬ 
er  frequencies.  The  optimum  ranges  of  diameters  for  stability,  however,  are 
approximately  the  same  in  both  cylindrical  and  one-dimensional  geometries. 

For  the  first  axial  and  tangential  modes,  however,  the  trends  are  significant¬ 
ly  different  from  those  of  the  axial  mode.  The  two-dimensional  cylindrical 
system  is  more  stable  than  the  one-dimensional  system  for  all  frequencies. 

There  is  an  indication  that  optimum  particle  diameters  shift  toward  larger 
sizes  as  the  frequencies  decrease.  For  further  details,  see  Reference  [13] 
or  Appendix  2. 
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to  the  gain  of  energy,  increasing  with  frequency. 

It  is  interesting  to  note  that  the  trend  of  convec¬ 
tion  into  domain  (E)  is  opposite  from  the  behavior 
of  the  surface  convection.  Note  that  energy  is 
lost  due  to  convection  into  the  domain  at  a  low 
frequency,  but  it  gains  at  a  higher  frequency,  con¬ 
trary  to  the  case  of  surface  convection.  This 
trend  appears  to  be  the  result  of  the  turbulent 
mean  flow  field. 

Acoustics-coupled  vorticity  instability  growth 
constants  (a  )  versus  vortical  frequencies  for  the 
laminar  flow,  with  ■  2553  Hz,  Re*l03,  6”  I  *♦ " ,  are 
shown  in  Fig.  5.  In  view  of  negligible  positive 
growth  constants,  it  is  concluded  that  the  system 
does  not  appear  to  be  unstable.  If  turbulent  flow 
is  considered,  however,  several  unstable  motions  at 
low  frequencies  are  observed  (Fig.  6).  As  the 
transition  angle  increases  to  8»20“,  the  system 
gradually  turns  to  instability  at  higher  frequenc¬ 
ies  for  laminar  flow  (Fig.  7).  Such  instability 
appears  to  occur  at  a  lower  frequency  if  turbulent 
flow  is  considered  (Fig.  8).  These  trends  are  more 
pronounced  as  the  transition  angle  is  increased  to 
0»36*.  To  compare  these  observations  with  the  case 
of  a  low  acoustic  frequency  combination,  the  result 
of  <i)N  «  36  Hz  is  investigated  (Fig.  II).  It  is 
clear  that,  for  a  small  transition  angle  ( 9->  1  *4 * )  , 
positive  growth  constants  increase  in  magnitude 
significantly,  not  only  at  low  frequencies,  but  al¬ 
so  at  high  frequencies  in  contrast  to  the  case  of 
high  acoustic  frequency.  »  2553  Hz,  as  shown  in 
Fig.  6. 

The  foregoing  discussions  on  acoustics-coupled 
vortical  instability  lead  us  to  a  critical  point  of 
re-examination  of  the  present  theory.  It  is  as¬ 
serted  that  there  exists  a  combination  of  any  a- 
coustic  frequency  with  any  other  vortical  frequency 
which  may  produce  a  certain  state  of  stable  or  un¬ 
stable  motion.  It  is  possible  that  not  all  combin¬ 
ations  of  acoustic  and  vortical  frequencies  theore¬ 
tically  postulated  may  be  excited.  In  fact,  only  a 
limited  number  of  combinations  would  be  considered 
significant  in  practice.  To  determine  whether  any 
of  the  combinations  of  acoustic  and  vortical  fre¬ 
quencies  are  excited,  one  may  resort  to  the  plots 
of  both  acoustic  and  vortical  modes  of  all  possible 
combi nat ions ,  the  process  of  which  can  easily  be 
automated  by  means  of  computer  graphics.  Such  an 
effort  is  currently  under  progress.  As  a  result  of 
this  analysis,  it  is  possible  to  construct  stabili¬ 
ty  boundaries  similar  to  those  via  solutions  of  the 
Or r-Sommerfe I d  equation.  It  is  seen  that  for  <aN  « 
2553  Hz,  e-Ht’,  the  stability  boundaries  for  acous- 
t ic-vort ic i ty  interactions  with  various  vortical 
frequencies  versus  Reynolds  numbers  exhibit  multi¬ 
ple  islands  as  shown  in  Fig.  12  for  the  laminar 
flow.  If  turbulent  flow  is  considered,  however, 
the  stability  boundaries  expand  significantly  In 
size,  and  move  toward  lower  vortical  frequencies. 

It  is  shown  that  the  critical  Reynolds  number  ap¬ 
pears  to  be  around  1*00-500.  The  most  interesting 
aspect  of  the  mu  1 1  i  p  i  e,-  i  s  I  and  stability  boundaries 
is  that  large  bay  areas,  which  indicate  stable  re¬ 
gions,  exist  for  certain  combinations  of  vortical 
frequencies  and  Reynolds  numbers.  The  trends  ob¬ 
served  herein  are  based  on  the  limited  number  of 
data  reductions.  It  is  anticipated  that  more  con¬ 
clusive  observations  will  be  made  available,  pend¬ 
ing  additional  computer  graphics  data  reductions. 


6.  Conclusions 

There  are  several  major  findings  in  this  work. 
They  are  summarized  as  follows: 

(1)  A  consistent  derivation  of  stability  inte¬ 
grals  for  the  growth  constants  associated  with 
acoustic  and  vortical  oscillations  leads  to  var¬ 
ious  terms  which  have  appeared  for  the  first 
time.  As  a  special  case,  the  flow  turning  term 
arises  as  a  consequence  of  integrations  by  parts 
twice  of  the  convective  terms. 

(2)  The  K-e  turbulence  model  appears  to  provide 
reasonable  flow  fields  consistent  with  the  tran¬ 
sition  angles  of  a  circular  cross-section,  which 
are  used  in  the  calculation  of  growth  constants. 

(3)  The  convection  into  domain  at  a  low  fre¬ 
quency  leads  to  the  loss  of  energy,  a  trend  re¬ 
versed  at  a  higher  frequency  resulting  in  the 
gain  of  energy.  This  phenomenon  Is  contrary  to 
the  case  of  surface  convection,  which  Is  consid¬ 
ered  to  be  the  effect  of  the  turbulent  mean  flow 
field. 

(4)  Unstable  motions  at  low  frequencies  appear 
'  to  be  the  result  of  turbulent  flow.  Increasing 

In  ma^iitude  for  larger  transition  angles. 

(5)  There  exists  a  combination  of  any  acoustic 
frequency  with  any  other  vortical  frequency 
which  may  produce  a  certain  state  of  stable  or 
unstable  motion.  Excitations  of  such  frequenc¬ 
ies  may  not  be  assured  for  all  combinations. 

(6)  Stability  boundaries  similar  to  those  as 
calculated  by  the  Orr-Sommerfeld  equation  are 
constructed  for  a  given  acoustic  frequency,  re¬ 
sulting  in  multiple  islands. 

(7)  The  effect  of  turbulent  flow  on  stability 
boundaries  is  to  expand  their  sizes  and  to  move 
toward  lower  vortical  frequencies,  with  large 
bay  areas  indicating  stable  regions  for  certain 
combinations  of  vortical  frequencies  and  Rey¬ 
nolds  numbers. 
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Simultaneous  iterative  solutions  of  Eqs. (67, 
71)  provide  the  turbulent  mean  flow  field.  To¬ 
gether  with  the  eigenvalue  analyses  for  the  acous¬ 
tic  field  of  Eq.  (53)  and  the  oscillatory  vortical 
field  given  by  Eq .  (58),  the  turbulent  mean  flow 
calculations  will  then  provide  necessary  informa¬ 
tion  for  the  calculation  of  growth  constants  by 
means  of  Eqs.  (33*35).  The  stability  integrals  as 
dictated  by  Eqs.  (3k-35)  can  be  performed  ideally 
by  the  Gaussian  quadrature  approximations  (251 • 


5.  Numerical  Applications 


formulations  presented  in  the  previous  sections,  a 
typical  solid  propellant  rocket  motor  with  an  axi- 
symmetric  cylindrical  geometry  has  been  analyzed. 
As  shown  in  Fig.  1,  k8  linear  isoparametric  ele¬ 
ments  with  63  global  nodes  are  chosen  to  model  the 
cavity  geometry.  Transition  angles  for  the  cross- 
section  of  the  motor  are  varied  G=  14°,  20°,  3k°. 

Various  constants  used  in  this  analysis  in¬ 
clude  the  following: 

Admittance,  AM  *  10.0  applied  normal  to  the 
burning  surface  for  the  region  designated  by  a 
through  c  (Fig.  1) 

Empirical  constants  in  turbulent  mean  flow 

C  =  0.09 
b 

a,  -  1.0 

k 

o  -  1.208 

V  ’-8 

C  -  0.6 

r2 

Reynolds  number,  Re  *  1 0 3 ,  10",  10s 

The  results  of  mean  flow  calculations  are 
shown  in  Figs.  2 -A.  In  general,  turbulent  velo¬ 
city  profiles  are  steeper  in  the  vicinity  of  the 
wall,  and  flatter  in  the  core  region  of  the  com¬ 
bustion  chamber  than  laminar  velocity  profiles. 
These  variations  are  caused  by  the  shear  stresses 
which  are  increased  by  the  turbulent  effects  (29], 
a  similar  trend  as  the  turbulent  pipe  flow  (30]. 
However,  it  should  be  noted  that  the  recirculation 
which  would  occur  near  downstream, i s  greatly  sup- 
ressed. since  the  admittance  is  applied  normal  to 
the  burning  surface,  contrary  to  the  usual  en¬ 
trance  boundary  condition  of  a  pipe  flow  in  which 
the  boundary  velocity  is  applied  parallel  to  the 
axis.  Separation  and  recirculation  flows  appear 
in  both  the  laminar  and  the  turbulent  cases  as  the 
transition  angles  are  increased.  But,  in  a  turbu- 
'ent  flow,  they  appear  only  at  a  relatively  large 
transition  angle  since  an  increase  of  shear 
stresses  with  turbulence  near  the  wall  overcomes 
the  effect  of  adverse  pressure  gradients. 

Typical  results  for  the  growth  constants, 
based  on  the  turbulent  mean  flow,  are  shown  in 
table  1,  listing  contributions  of  various  terms  in 
the  stability  integrals  for  the  acoustic  frequen¬ 
cies  (<uu  =  3k  Hz,  k76  Hz)  in  combination  with  the 

vortical  frequency  «  2k  Hz  for  Reynolds  number 
Re»l03,  and  transition  angle  0”lk°.  Note  that  the 
so-called  pressure  coupling  and  velocity  coupling 
defined  as  the  surface  combustion,  designated  as 
(A),  provide  a  significant  source  of  acoustic  in¬ 
stability  (a^  ■  I5.A58)  for  *  3k  Hz,  whereas 

for  a  higher  acoustic  frequency  (u;  =  k 76  Hz),  the 

acoustic  growth  constant  assumes  *  -3k. 255.  an 
indication  of  strong  stability.  The  flow  turning 
effect  which  is  included  in  the  surface  convec¬ 
tion,  designated  as  (B)  ,  exhibits  a  similar  trend 
as  the  surface  combustion,  but  less  in  magnitude. 
For  the  Reynolds  number  and  the  transition  angle 
investigated  (Re=IOJ,  =lk°),  the  effect  of  vis¬ 
cosity  is  negligible  as  indicated  by  (C),  (F) ,  and 
(G) .  However,  it  will  be  found  in  the  later  anal¬ 
ysis  that  there  is  a  significant  contribution 
from  the  viscous  terms  through  the  mechanism  of  a 
vorticity  generation  in  creating  the  Reynolds  num¬ 
ber-dependent  stability  boundaries.  The  combus¬ 
tion  into  the  domain,  as  identified  by  (0),  leads 


To  verify  the  validity  of  the  theory  and  the 


(61) 


tions  2  and  3. 

k.  Finite  Element  Analysis 


The  use  of  finite  elements  in  fluid  mechanics 
problems  has  increased  significantly  in  recent 
years  [25).  Especially,  finite  elements  offer 
greater  flexibilities  for  the  complex  geometry 
such  as  combustion  chambers  in  the  solid  propel¬ 
lant  rocket  motors.  In  this  paper,  finite  ele¬ 
ments  are  used  in  three  parts,  i.e.,  the  eigenval¬ 
ue  analysis  in  the  classical  acoustics  and  the  os¬ 
cillatory  vortical  field,  calculations  of  mean  ve¬ 
locity  fields,  and  stability  integrals. 

first  of  all,  we  return  to  a  classical  acous¬ 
tic  problem  characterized  by  Eqs.  (2k-25)  .  Using 
Galerkin  finite  elements.  Eg.  (2k)  leads  to 


<Vii  +kNPN>Vn 


(53) 


where  #  is  the  test  function  which  is  set  equal 


01 


to  the  trial  function  such  that 


VV  "  WPNa 


(5k) 


Eq.  (53)  can  be  solved  by  the  finite  element  ei¬ 
genvalue  equation  of  the  form 


| A  .  -  k*B  I 
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where  k^  is  the  wave  number  from  which  the  acous¬ 
tic  frequency  c>N  may  be  determined. 


a8 


'id 


f  <t>  ,4>  . 

n  t  * 1  6 , 1 

f  $ 

in  01 


dfl 


(56) 

(57) 


Similarly,  the  vortical  fluctuation  Eq.  (51)  is 
cast  in  the  form 
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where  k  denotes  the  vortical  wave  number  from 


which  the  vortical  frequency  oil  may  be  calculated. 
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As  a  result,  frequencies  of  normal  acoustic 
modes  and  vortical  oscillations  can  be  obtained  by 
taking  the  real  parts  of  the  eigenvalues  of  Eqs. 
(55,58),  respectively.  The  vorticity  growth  con¬ 
stants  can  be  obtained  by  taking  the  imaginary 
parts  of  the  eigenvalues  of  Eq.  (58)  [26-291.  How¬ 
ever,  the  imaginary  parts  of  the  eigenvalues  of  Eq. 
(55)  are  absent  because  the  acoustic  growth  constants 
do  not  exist  in  the  normal  modes. 

The  mean  flow  calculations  are  performed  under 
the  assumption  of  a  steady  state,  isothermal,  in¬ 
compressible  flow.  It  follows,  then,  that  the  gov¬ 
erning  equations  of  the  mean  flow  are  simplified  as 
foi lows : 
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where 


1 


1  ,  I 


Re  rr  Re  Re. 
eff  t 


(6k) 


With  an  assumption  of  incompressible  flow,  the 
K-e  equations  can  be  simplified  as  follows: 
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To  solve  these  non-linear  equations  by  Galer¬ 
kin  finite  elements,  the  Newton-Raphson  method  is 
used.  In  order  *o  reduce  computational  efforts, 
the  Jacobian  matrices  for  the  flow  field  and  the 
K-e  system  are  separated,  which  will  then  be  updat¬ 
ed  iteratively  between  the  two  Jacobian  matrices. 
The  Newton-Raphson  scheme  for  the  velocity  field  is 
given  by 


j  (n)  ,x (n+i )  _  (n) 

pq  q  =  p 


(67) 


where 


X (n+1 >  ,  x(n)  *  AX(n+,) 


(68) 


with 


J(n)  -  f  0>  .6..^  0  .  +  *  < 

pq  lrt  a  Bj  *  k  y  >  j  a 

.  I 


'T  ,  u 
l.k  Yi  b 


X 

+  — 


1  *i,i.;,k  Reeff  i.jb.j  ik 


t  .i.Jd"  (69) 


.(n) 


f  ( t  T ,  .  a . . .  U  .  +  -  :  .  u  . . 

f  .1.  .u..)d'  -  f  (  —  u .  .n.i 
' .  J  ■  .  J  -  1  J-  '■  '  .J  1  > 

u.  .n.i-  )rt:'  (70) 


I 


Re 


eff 

I 


Re 


eff 


14 


aCu>  P  -  P/P  ,  u  »  u/a, 

rf  O 


(yP0/p) 


Substituting  these  relations  and  Eq.  (31)  to  Eq. 
(3A),  collecting  the  first  two  terms  of  (A)  and  the 
flow  turning  term,  and  rewrit  ing  in  the  same  form 
as  Eq.  (27),  we  obtain 

(kJ  -  k2)E2  =  ioakN  fr<u'PN  +  u)-ndf 


j  (57PN)20u*ndr 


It  is  seen  that  pu*n  is  equivalent  to  in  terms 

of  the  notations  used  in  [A],  Notice  that  Eq.  (36) 
above  is  the  same  as  Eq.  (A.1A)  in  [A]  neglecting 
the  particle  distribution  of  the  two  phase  flow. 


3.  Turbulent  Mean  Flow  And 
Fluctuation  Vorticities 

A  glance  at  Eq.  (33)  with  details  given  by 
Eqs.  (3A-35)  indicates  that  the  flow  field,  in¬ 
cluding  the  mean  velocity  and  the  vorticity  to¬ 
gether  with  the  vortical  component  of  the  fluctu¬ 
ation  velocity,  as  well  as  the  acoustic  pressure 
modes,  must  be  calculated.  To  this  end,  we  include 
turbulent  effects  on  the  mean  flow  field,  but  at 
this  time,  shock  waves  are  excluded  from  considera¬ 
tion.  Turbulent  modelling  for  solid  propellant  rock¬ 
et  motors  has  been  the  subject  of  controversy  parti¬ 
cularly  due  to  complex  burning  surface  phenomena. 
Pending  development  of  an  adequate  model  in  the  fut¬ 
ure,  we  examine  here  the  K-E  model  for  a  computatio¬ 
nal  purpose.  Thus,  the  Reynolds  decompos i t ion  becomes 

u  i  ( x  j ,  t )  2  u .  (x  j )  +  uTUj.t)  (37) 

Then,  the  time-averaged  velocity  is  defined  as 
u.  «  ouj/o  (38) 


(cu.)  =  0  (39) 

1 

It  follows  that 

pu.u.  *  pu.u.  +  (pu.)'u."  (AO) 

i  J  i  J  1  J 

where  the  bars,  tildes  and  primes  denote  time-aver¬ 
aged,  mean,  and  fluctuating  values,  respectively. 
Eqs.  (37-AO)  refer  to  the  time-averaged  and  fluc¬ 
tuating  components  of  the  scalar  fields,  P  and  T. 

Substituting  Eqs.  (37-AO)  into  Eqs.  ( I  -  A) , 
the  time-averaged  governing  equations  are  written 
as  fol lows : 
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In  terms  of  the  eddy  viscosity  hypothesis  (30],  the 


Reynolds  stress  -  (ou. ) 'u 7  can  be  expressed  as 

’  (puiruj  =  fk~t  Si.j  ( 
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tot 

with  p(  being  the  eddy  viscosity.  The  eddy  viscos¬ 
ity  can  be  expressed  as  a  function  of  turbulent 
kinetic  energy  K  and  energy  dissipation  E  by  the 
Prandt I -Kolmogorov  law 

tb7  -  V  r  (,,7) 

where  C  is  the  empirical  constant. 

P 

The  governing  equations  for  K  and  e  are  deriv¬ 
ed  from  the  momentum  equation  [  2A]  as 
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3t  1  ,1  Retak  ,J  ,J 
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in  which  the  following  non-d imens iona 1  quantities 
are  used: 

K  -  K/a2  ,  E  »  "L/a3,  Re£  -  PQaL/ut 

Furthermore,  and  o£  are  the  effective  Prandtl- 

Schmidt  numbers,  and  Cei  and  Cj.  are  the  empirical 
constants. 

In  addition  to  the  turbulent  mean  flow,  we  re¬ 
quire  that  the  vortical  component  of  the  fluctua¬ 
tion  velocity  and  the  fluctuation  vorticity  be  cal¬ 
culated.  In  this  regard,  we  set  u  .  «  XT.  +  £u*. 

P"l,  and  p“!  to  obtain  the  fir-t  order  perturba¬ 
tion  momentum  equation  in  the  form 
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It  follows  from  Eq .  (17)  and  Eq.  (50)  that 
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The  fluctuation  vorticity  is  then  calculated 


u  52 

1  1  jK  k  ,j 

in  which  the  vortical  component  of  the  fluctuation 
velocity  u  is  obtained  as  a  consequence  of  Eq. 
(51).  klJ 

In  the  sequel,  we  discuss  finite  element  so¬ 
lutions  to  all  of  the  equations  presented  in  Sec- 
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where  the  superscripts  (R)  and  (i)  refer  to  the 
real  and  the  imaginary  parts,  respectively. 

The  normal  velocity  at  the  surface  can  be  ex¬ 
pressed  in  terms  of  the  admittance,  A  -  A^  + 
and  the  mean  flow  Mach  number  R  such  that 

u.'n.  *  (U.'(R)  +  iu.'^Jn.  -  MiA P  /y  (31) 

III  I  I  N 

whereas  the  acoustic  fluctuation  velocity  in  the 
domain  is  given  by 

u.'  -  ~  Pu  .  (32) 

i  kN>  N ,  i 

For  the  purpose  of  investigating  the  coupling 
mechanism  of  acoustic  and  hydrodynamic  instabili¬ 
ties,  it  is  convenient  to  separate  Eq.  (30)  into 
two  parts : 

‘  ~  *  C  (33) 

where  and  refer  to  the  acoustic  and  the  hy¬ 
drodynamic  contributions,  respectively. 
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Note  that  the  first  term  of  the  boundary  in¬ 
tegrals,  designated  as  (8) 
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is  seen  to  be  identical  to  the  three-dimensional 
equivalent  of  the  one-d i mens iona 1  flow  turning 
term  which  appeared  in  Culick  [31,  but  which  did 
not  arise  in  due  course  of  mathematical  deriva¬ 
tions  for  the  three-dimensional  case  [Culick,  k). 
To  reconstruct  the  present  results  in  terms  of  the 
same  dimensions  and  notations  as  in  [k],  we  set 
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The  boundary  condition  can  be  obtained  from  the  mo¬ 
mentum  equation  by  constructing  a  dot  product  with 
the  normal  vector, 
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The  oscillatory  motions  in  both  the  acoustic 
and  the  vortical  fields  are  modeled  by 
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where  k  is  the  complex  dimensionless 
en  by 

frequency  giv 

k  *  w-ia 

(18) 

Here,  the  imaginary  part  a  is  known  as  the  growth 
rate. 

Substituting  Eqs.  05-18)  into  the  wave  equa¬ 
tion,  we  arrive  at  the  nonhomogeneous  Helmholtz  e- 
quat ion , 

P  . .  +  k2P  -  h 
»  1  1 

(19) 

and  the  boundary  condition, 
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-p  .n.  -  f. 
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The  nonhomogeneous  terms  h  and  f  are  given  by 
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and  Eu  is  given  by 

N 
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CN 
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(26) 

At  this  point,  an  important  remark  is  in  or¬ 
der.  A  close  examination  of  Eq.  (23)  reveals  that 
the  domain  integral  on  the  right  hand  side  of  Eq. 
(23)  contains  the  terms  from  the  momentum  equation 
which  were  differentiated  once.  Thus,  these  terms 
must  be  integrated  by  parts  to  produce  an  acoustic 
boundary  condition.  The  resulting  domain  integral 
represents  the  functional  space  equivalent  to  the 
mean  flow  characterized  by  the  Navier-Stokes  sys¬ 
tem.  This  implies  that  an  additional  integration 
by  parts  must  be  carried  out  such  that  the  familiar 
Neumann  boundary  conditions  may  be  brought  to  the 
surface.  The  boundary  integrals  arrived  at  in  this 
manner  account  for  the  stress  and/or  the  pressure 
boundary  conditions  by  means  of  velocity  gradients. 
Such  Neumann  boundary  conditions  stem  from  the  con¬ 
vective  and  viscous  terms. 

In  view  of  these  requirements,  the  resulting 
expression  upon  integration  by  parts  twice  of  Eq. 
(23)  takes  the  form 
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that 


Squaring  both  sides  of  Eq.  (18)  and  noting 


Making  use  of  the  Green's  function  Integral 
(23).  it  can  be  shown  that 

(kJ  *  kN)EN  ■  l  "PNd“  +  Jp  EVr  (2*> 

where  the  unperturbed  mode  shape  P,,  and  its  fre- 
quency  k„  are  determined  from  the  classical  acous- 
tic  problem: 
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Since  t‘  <<  [ 2 ,  from  the  condition  set  by  Eq. 
(28),  it  is  now  possible  to  solve  for  the  growth 
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stability  can  be  significant  [7-13].  Although  it 
can  be  argued  that  the  hydrodynamic  instability  may 
not  occur  in  high  Reynolds  numbers,  the  turbulent 
shear  layer  instabilities  have  been  found  to  be  af¬ 
fected  by  various  combinations  of  Strouhal  numbers 
and  Reynolds  numbers.  The  acoustic  field  may  in¬ 
teract  with  vortex  motions  known  as  the  "feedback" 
resulting  in  the  votex  generated  sound  flA],  Some 
studies  [15,16]  indicate  that  the  vortices  may  un¬ 
dergo  "clipping",  a  phenomenon  corresponding  to  the 
vortex  disruption.  It  is  also  possible  that  later¬ 
al  periodic  motion  of  the  vortex  street  known  as 
"jitters"  may  lead  to  partial  or  complete  escape  of 
the  vortices  [15,16].  Whether  these  conditions 
prevail  in  large  rocket  motors  in  which  flow  sepa¬ 
rations  at  interface  restrictors  or  inhibitors  are 
likely  to  produce  vortex  motion  must  be  clarified. 
No  simple  models  such  as  hyperbolic  tangent  velo¬ 
city  profile  for  the  shear  layer  [11,17]  and  tem¬ 
poral  or  spatial  growth  theories  [18,19]  appear  to 
be  adequate  for  the  interactions  of  acoustic  and 
vortical  oscillations  in  a  rocket  chamber. 

In  the  previous  papers  [20-22],  finite  element 
applications  to  the  combustion  instability  analysis 
were  discussed.  Although  a  rigorous  mathematical 
formulation  of  the  stability  integral  was  present¬ 
ed,  the  mean  flow  calculations  did  not  include  tur¬ 
bulent  flows.  Since  the  turbulent  flow  field  is 
involved  in  shear  boundary  layers  and  vortex  mo¬ 
tions,  it  is  intended  that  this  subject  be  consid¬ 
ered  in  the  mean  flow  calculations  and  subsequent¬ 
ly  in  the  stability  integral.  Shock  waves  will  not 
be  included  in  the  present  paper. 

The  numerical  results  for  certain  combinations 
of  acoustic  and  vortical  frequencies  indicate  that 
stability  boundaries  for  acoust i cs-coupl ed  hydro¬ 
dynamic  oscillations  are  somewhat  similar  to  the 
classical  hydrodynamic  stability  boundaries,  but 
they  occur  in  the  form  of  multiple  islands.  The 
turbulent  flow  field  appears  to  contribute  toward 
instability,  and  this  trend  increases  with  larger 
transition  angles  of  the  rocket  motor  cross-sec¬ 
tion. 

In  what  follows,  pertinent  governing  equations 
are  presented,  from  which  the  expression  for  the 
growth  constant  coupling  the  acoustic  and  hydrody¬ 
namic  oscillations  is  derived  in  section  2.  Subse¬ 
quently,  the  K-e  turbulence  model  and  calculations 
of  the  vortical  component  of  fluctuation  velocities 
and  the  fluctuation  vorticity  are  described  in 
section  3.  It  is  also  shown  in  section  A  that  the 
finite  element  is  one  of  the  most  expeditious  me¬ 
thods  of  calculation.  Numerical  results  via  finite 
elements  are  then  displayed  and  pertinent  discus¬ 
sions  and  conclusions  are  presented  in  sections  5 
and  6,  respectively. 
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where  the  commas  denote  partial  derivatives,  the 
repeated  indices  imply  summing.  The  following  non- 
dimensional  quantities  are  used  in  the  above  equa¬ 
tions  : 

u.  *  u./a,  a  *  (yP  /o  ,  P  =*  P/p 

II  o  o  o 

T  -  cp(y-l)T/aJ,  x.  •  x./l,  t  -  at/L 

Re  »  0  aL/u  ,  p  -  p/p 
o  o 

where  the  double  bars  denote  dimensional  quanti¬ 
ties. 

Interactions  between  acoustic  and  vortical  os¬ 
cillations  can  be  introduced  by  superimposing  the 
acoustic  component  upon  the  vortical  component  of 
the  perturbed  velocity  in  the  form 

ui  -  ut  +  e<uf  ♦  u*)  (5) 

where  the  bars,  primes,  and  asterisks  indicate  the 
mean  flow,  the  acoustic  and  vortical  oscillations, 
respectively;  and  e  represents  the  perturbation 
parameter.  On  the  other  hand,  the  pressure  and  the 
density  are  given  by 
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the  vorticity  field  is  given  by 
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where  is  the  permutation  symbol.  In  view  of 

Eqs.  (1-10) ,  it  can  be  shown  that  the  nonhomogen- 
eous  wave  equation  takes  the  following  form: 


(ID 


The  basic  governing  equations  for  compressible 
viscous  flow  without  particle  distributions  are  re¬ 
presented  as  follows: 
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APPENDIX  1 


INFRACTIONS  OF  UNSTEADY  ACOUSTIC  AND  VORTICAL 
OSCILLATIONS  IN  AX  I  SYMMETRIC  CYLINDRICAL  CAVITY 

T.J.  Chung*  and  J.L.  Sohn** 
Department  of  Mechanical  Engineering 
The  University  of  Alabama  in  Huntsville 
Huntsvi I le,  Alabama 


Abstract 

A  physical  phenomenon  of  interactions  between 
the  acoustic  and  vortical  oscillations  is  examined 
herein.  This  subject  is  important  in  rocket  motor 
chambers  when  the  vorticity  field  is  coupled  with 
acoustic  pressure  oscillations.  In  the  past,  the 
acoustic  combustion  instability  was  studied  inde¬ 
pendently  of  the  hydrodynamic  instability  induced 
by  vortex  motions  and  turbulent  shear  boundary  lay¬ 
ers.  However,  it  is  quite  conceivable  that  these 
two  distinctly  different  oscillations  are  coupl¬ 
ed  and  interact  together  in  the  flow  field  of  a  so¬ 
lid  propellant  rocket  motor.  The  present  paper  in¬ 
troduces  an  analytical  approach  to  resolving  the 
seemingly  complex  phenomena  of  mutual  interactions 
between  the  acoustic  and  vortical  oscillatory  mo¬ 
tions.  Toward  this  end,  governing  equations  for 
all  variables  are  constructed,  and  finite  elements 
are  applied  to  solye  the  governing  equations.  Com¬ 
bustion  instability  integrals  including  the  mean 
flow  field,  perturbed  acoustic  oscillations,  and 
oscillatory  velocities  and  vortices  are  also  deriv¬ 
ed  and  calculated  by  finite  elements.  From  the 
growth  constants  for  acoustic  and  hydrodynamic  con¬ 
tributions,  stability  boundaries  are  determined  in 
terms  of  Reynolds  numbers.  The  numerical  results 
indicate  that  an  overall  Instability  phenomenon  re¬ 
sults  from  certain  combinations  of  acoustic  and 
vortical  frequencies.  It  is  also  found  that  stabi¬ 
lity  boundaries  for  acoustics-coupled  hydrodynamic 
oscillations  are  somewhat  similar  to  the  classical 
hydrodynamic  stability  boundaries,  but  they  occur 
in  the  form  of  multiple  islands.  The  turbulent 
flow  field  appears  to  contribute  toward  instability 
and  this  trend  increases  with  larger  transition  an¬ 
gles  of  the  rocket  motor  cross-section. 

Nomenclature 

A  admittance  at  the  burning  surface 

a  sonic  velocity  (  Ypq/do  ) 

c  specific  heat  at  constant  pressure 

P 

t  length  of  the  combustion  chamber 

k  complex  wave  number  (w-ia) 

K  turbulence  kinetic  energy 

M  Mach  number 

P  pressure 

Re  Reynolds  number  (poaL/u) 

T  temperature 

t  time 

u .  velocities 

x.  spatial  coordinates 

i 
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a  growth  constant 

£  turbulent  energy  dissipation, 

perturbation  parameter 
e ; j k  permutation  symbol 

a  effective  Prandtl-Schmidt  number 

$3  finite  element  interpolation  function 

vorticities 
p  density 

Y  specific  heat  ratio 

y  viscosity 

T  boundary 

ft  domain 


Subscripts  and  Superscripts 


acoustic  fluctuation 
vortical  fluctuation 
turbulence 
norma  1  mode 
acoustic  field 
hydrodynamic  field 


I.  Introduction 

The  flow  field,  such  as  occurs  in  solid  pro¬ 
pellant  rocket  motors,  offers  a  fertile  ground  for 
fundamental  research  in  fluid  mechanics  and  heat 
transfer.  Combustion  induces  not  only  the  mean 
flow  field,  but  also  acoustic  pressure  oscillations 
and  possibly  vortex  fluctuations  together  with  tur¬ 
bulent  shear  boundary  layers.  Furthermore,  shock 
waves  are  commonplace  in  most  instances.  Obvious¬ 
ly,  a  most  rigorous  analysis  taking  into  account 
all  of  these  phenomena  would  be  difficult,  if  not 
impossible.  However,  with  the  advent  of  the  elec¬ 
tronic  computer  and  the  modern  technology  of  numer¬ 
ical  methods,  it  has  become  feasible  to  resolve 
hitherto  unsolved  problems. 

Despite  difficulties  in  analytical  and  numeri¬ 
cal  solutions  to  the  complex  physical  phenomena  in 
a  rocket  motor  chamber,  many  researchers  have  con¬ 
tributed  to  the  advancement  of  analysis  and  design 
of  successful  rocket  motors.  A  large  body  of  lit¬ 
erature  exists  relative  to  this  subject,  the  study 
of  which  has  been  pioneered  by  Crocco  [I],  Cantrell 
and  Hart  [2],  Culick  [3,4],  and  others.  Flandro 
and  Jacobs  [5],  among  others,  have  noted  that  vor¬ 
tex  shedding  may  lead  to  an  instability  in  solid 
propellant  rocket  motors.  It  is  quite  possible 
that  high  speed  mean  flows  also  affect  the  stabili¬ 
ty  [6]  significantly. 

The  basic  mathematical  formulations  of  combus¬ 
tion  instability  were  contributed  by  Culick  [3,4], 
Recently,  it  has  been  observed  in  both  full-scale 
firings  and  cold  flow  simulations  that  interac¬ 
tions  of  acoustic  and  hydrodynami ca 1  (vortical)  in- 
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SECTION  4 


UNSTEADY  RESPONSE  OF  BURNING  SURFACE  IN 
SOLID  PROPELLANT  COMBUSTION 

4 . 1  Summary 

The  response  function  as  related  to  the  combustion  of  solid  propellants 
has  been  studied  by  numerous  investigators.  Some  of  the  review  articles  in¬ 
clude  Culick  [1]  on  homogeneous  propellants  [2-6]  and  Cohen  [8]  on  heterogen¬ 
eous  propellants  [9-14] .  Both  of  them  are  limited  to  one-dimensional  and 
quasi-steady  situations.  In  recent  years,  some  works  toward  non-steady  prob¬ 
lems  have  been  attempted  [15-17].  The  majority  of  the  discussions  of  response 
functions  in  the  literature  are  concerned  with  pressure-coupling  usually  appli¬ 
cable  in  the  linear  stability.  Computations  of  the  response  function  for  velo¬ 
city-coupling  important  in  a  nonlirtear  process,  however,  remain  in  a  state  of 
infancy,  although  some  initial  attempts  toward  this  subject  have  been  made  [16, 
18]. 

As  noted  by  Flandro  [16],  the  velocity-coupling  may  be  accommodated  by 
second  order  perturbations  of  multi-dimensional  gas /solid  governing  equations. 
This  approach  deviates  drastically  from  the  one -dimensional  analysis  which  has 
been  adopted  for  nearly  three  decades.  With  modem  computers  and  various  tools 
of  numerical  analysis  available,  however,  it  seems  possible  to  relax  many  un¬ 
desirable  restrictions.  Observations  indicate  that  combustion  oscillations 
are  time-dependent  and  often  nonlinear  as  influenced  by  turbulent  flowfields, 
which  may  lead  to  erosive  burning  and  unstable  oscillations.  Unfortunately, 
however,  a  computational  tool  for  the  most  exact  analysis  of  complicated  phys¬ 
ical  phenomena  such  as  these,  even  if  developed,  will  not  fit  into  the  current¬ 
ly  available  computer.  Thus,  naturally,  some  approximations  and  simplifica¬ 
tions  must  still  be  introduced  to  any  theoretical  formulations  conducive  to 
numerical  analysis. 

With  this  in  mind,  some  additional  materials  and  reassessments  concern¬ 
ing  the  multi-dimensional  calculations  of  combustion  response  functions  are 
presented  herein  as  an  extension  to  the  previous  paper  [17].  The  influence 
of  heterogeneity,  surface  roughness,  particulate  matter,  or  turbulence  will 
not  be  considered  at  this  time.  We  include  the  effect  of  radiation,  adopt  a 
simple,  premixed,  single-step  laminar  flame,  expand  the  solid-gas  governing 
equations  into  first  and  second  order  perturbations,  and  finally  perform  ei¬ 
genvalue  analyses.  Complicated  boundary  conditions  at  the  solid-gas  inter¬ 
face  and  flame  edges  are  imposed  ideally  by  means  of  Lagrange  multipliers 
together  with  finite  elements.  Calculations  are  carried  out  for  various  in¬ 
cident  angles  of  impressed  pressure  waves  at  the  flame  edge  boundaries  [16]. 

Computed  results  show  that  the  natural  frequencies  are  clustered  around 
low  frequency  range  (ca<10) .  Spatial  distributions  of  field  variables  cor¬ 
responding  to  the  computed  frequencies  indicate  that  the  oscillatory  behav¬ 
ior  is  pronounced  at  upstream  and  gradually  diminishes  toward  downstream.  Re¬ 
sponse  functions  oscillate  in  the  axial  direction  with  peaks  occurring  at  the 
midstream,  but  diminishing  toward  flame  edges  for  both  first  and  second  or¬ 
der  perturbations.  It  is  also  shown  that  two-dimensional  response  functions 
are  multi-peaked  and  may  become  negative  as  energy  sinks.  Finally,  it  is 
seen  that  the  effect  of  radiation  is  more  pronounced  in  the  second  order  per- 
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Abstract  j 

The  paper  discusses  the  formulation  and  solu-  ! 
tion  strategies  to  analyze  combustion  instability 
due  to  aluminum  particle  distributions  in  the  solid 
propellant  rocket  motor  chambers.  Specifically,  the 
finite  element  method  is  utilized  in  order  to  ac¬ 
commodate  complicated  geometries  and  boundary  condi¬ 
tions.  To  demonstrate  the  validity  ,  one  dimension¬ 
al  results  are  first  compared  with  those  of  analyti-j 
cal  solutions.  A  similar  process  is  then  extended  t 
to  handle  two-dimensional  problems. 
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Nomenclature 


I 


speed  of  sound  | 

specific  heat  of  particle  ! 

specific  heat  of  gas  at  constant  pressure,! 
volume  ! 

specific  heat  of  mixture  at  constant  pres-' 
sure,  volume 

particle  diameter  ' 

total  energy  of  gas 
total  energy  of  particle 

drag  force  j 

wave  number 

vector  normal  to  surface 
mean  pressure 
acoustic  pressure 

Prandtl  number  ■ 

gas  constant 

temperature  of  gas 

temperature  of  particle 

radial  coordinate 

gas  velocicy 

particle  velocity 

particle  mass  rate  of  flow  i 

axial  coordinate  I 

growth  rate 

ratio  of  particle  density  to  gas  density 
specific  heat  ratio  of  gas 
specific  heat  ratio  of  mixture 

r  I 

total  density 
density  of  gas,  particle 
tangential  coordinate 
ratio  R/Cv 

frequency  ! 

thermal  conductivity 
viscosity  constant 

finite  element  interpolation  function 
vorticity  vector 
dynamic  relaxation  time 
thermal  relaxation  time 


1.  Introduction 

The  study  of  acoustic  energy  losses  due  to 
aluminum  particles  in  the  solid  propellant  rocket 
motor  combustion  chamber  has  been  carried  out  by 
a  number  of  investigators.  Epstein  and  Charhart  [1] 
studied  the  absorption  of  sound  in  suspensions  of 
non-interacting  inert  spherical  particles  and  uni¬ 
form  temperatures.  The  validity  of  this  investiga¬ 
tion  was  subsequently  substantiated  by  other  re¬ 
searchers  [  2-4  )  .  In  a  rocket  motor,  the  particles 
are  neither  spherical  nor  inert  and  subjected  to 
nonuniform  temperature  distributions  [5). 

! 

Despite  the  extensive  research  on  the  subject 
of  acoustic  energy  dissipation  due  to  particle  damp¬ 
ing  in  the  rocket  motor  [5-10],  calculations  of  the 
stability  integral  arising  from  particle  damping 
are  limited  to  simple  one-dimensional  cases. 

Thus,  the  purpose  of  the  present  study  is  to 
demonstrate  the  feasibility  of  mathematical  formula¬ 
tions  and  numerical  calculations  via  finite  elements 
Furthermore,  interactions  of  particle  damping  with 
fluid  viscosity  and  heat  transfer  are  included.  It 
is  shown  that  additional  boundary  and  domain  terms 
arise  from  integrating  by  parts  "twice"  of  the 
Green  function  stability  integral  containing  the 
momentum  equation.  A  simple  example  problem  is 
solved  for  comparison  with  the  analytical  solution, 
i 

2.  Coverning  Equations 

Conservation  li  s  for  the  mass,  momentum,  and 
energy  in  the  two-phase  flow  of  the  gas/particle  mix¬ 
ture  are  written  as  follows  [5]: 


Continuity 
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—  (pRc  +  Ppop)  +  7  •  t.yu  +  PpvV 


(u  )  -  Q  »  f- 
-P 


(4) 


where 


where  Wp  is  the  particle  mass  rate  of  flow  and 


S  =  p[V2u  +  -J  7  (7  -u)] 


Q  =  V  •  q^  +  7  ‘  q^  (6) 

2  ,  3ui  3ui  3u;  3u< 

(f  (7) 

(c)  (R) 

with  q  and  q  being  the  conduction  heat  flux 
and  radiation  fieat  flux,  respectively. 

Combining  (1),  (2)  and  (3),  we  obtain 


p  TjT  +  pg(u  -7)u  +  Vp  =  Fp  +  g  + 


!p  =  "  Ppt(lT  }  +  (“P  '  V)“p]  <9) 

C  *  (Up  -  u)wp  (10) 

Combining  (1),  (2),  and  (4)  leads  to 

cgCv  3t'  +  cgCv(“  '7)T  +  pV'“  ’  Q  +  ^  +  Qp 
-  u  •  S  +  (ep  -  e)wp  +  (up  -  u)Fp  +  u  •  O  (11) 


QP  ■  -pr^  +  %  '  7>tp 


Now  let  us  denote  that 


u  =  u  +  5u 
~P  -  -P 


T  =  T  + 

P  P 

In  view  of  (13)  we  rewrite  (9)  in  the  form 


pp[  -  +  (»  -v)  ui  +  «Fp 


r'F  =  -p  [  — +  (u  •  V)  ru  (■  (  u  •  V)u  (16) 
~P  p  3t  -  ,p  -P 


Substituting  (15)  into  (8)  yields 


+  f.(u  -y)u  +  7P  -  H  +  S 


6Qp  =  -PpC[-j^+  (u  •V)6Tp  +  (5u  -7)T  ]  (19) 


with  C  being  the  specific  heat  of  the  particle. 
Substituting  (18)  into  (11)  gives 


pg(Cv  +  Be)  +  Pg  (Cv  +  gc)  (u  -7>T  +  p7-  u 

=  Q  +  <5  +  Y  (20 

where  Cv  is  the  specific  heat  of  gas  at  constant 
volume  and 

Y  =  6Q  +  (e  -  e)w  +  (u  -  u)T  -  u  •  a 
P  p  p  ~p  ..  -p 


To  obtain  the  wave  equation  we  make  use  of 
the  perfect  gas  law 

p  *  p  RT  (22 

which  is  substituted  into  (20).  We  then  differ¬ 
entiate  (20)  with  respect  to  time  and  substitute 
from  the  spatial  derivative  of  (17).  Denote  the 
specific  heats  of  the  mixture  at  constant  volume 
and  constant  pressure,  respectively,  as 

Cv  +  3C  C  +  gC 


it  then  follows  that  the  gas  constant  for  the  mix¬ 
ture  takes  the  form 


R  -  C  -  C  =  -~-r 

P  v  1  +  g 


Substituting  (22)  and  (1)  into  (20)  yields 

+  (l+()pV  ’  u  +  (u  -7)p  =  £(Q  +  $  +  Z)  (23) 


where  c  =  R/C 

^  \7 


Z  «  Y  +  (1+B)C  Tw 
v  P 


New  the  time  derivative  of  (23)  and  spatial  deri¬ 
vative  of  (17)  are  written,  respectively,  as  fol¬ 
lows  : 


77?  +  (l+c)  [  •  U  +  p  ~  (V-  u)  1 

+  37  -  ",  37  (Q  +  *  +  z) 


8  P  g 


(1  f  ,-) »  B  *  o  /p 


P  (V  *  11 )  +  >3(u  *V)u  +  V?P 


=  V  •  H  +  V  •  S  +  B 


H  -  rF  +  ,"T 
P 

Similarly 

On  1  ‘PC  (  ~  +  (u  •  v)T  1  +  :Q 


(-  V  P  +  H  +  S) 


Substituting  (25)  into  (24)  gives 


and 


Z*  +  r>  Il>. 

3tfc  *  3t 


u  (  -p/  *(u*V)u  -  V'p 


+  7  •  H  +  7-S  +  BJ  +  —  (u •  7)p  -C  ^(Q+^+Z)  (26) 


For  linear  stability,  we  introduce  the  per¬ 
turbation  expansion  of  0(")  as  follows: 


u  *  u  +  i.  (u*  +  u*) 
P  =  P  +  "p’ 

c  =  c  +  cp' 


e  =  e  +  e1  ,w  -  w  +  ew 1  ,  Q  =  Q  +  ^0  * 
P  P  P 


1  =  1  +  {.(  t*  +  i  *) 

H  =  H  +  c(H'  +  H*) 

”  —  7  +  .;(>*'  +  a*) 


Z  +  EZ' 


where  the  bar,  astrisk,  and  prime  denote  the  mean 
flow,  acoustic  fluctuation,  and  vortical  fluctua¬ 
tion,  respectively. 


Now,  collecting  the  terms  of  0(r)  ,  we  have 


r!  n'  -L  92p'  _  h 


3t  ‘ 


(27) 


where 


h  =  -  n  (72(u-u')  +  V2(u-u*)  -  v  •( u  x  r 

+  u*  x  |  +  u'  x  c  )  1  +  f2lf  7-u 


r* 


+  4-  u-  (Vp)  +  7-H'  +  7-S  ' 

a  „  r)  t  — 


a2 

9 


a1  3t 


(Q’+i’+Z') 


(28) 


where 


-,  --^  -  Po  y  Po 

a  -  vRT  =  V  -z-  =  7%  ,  - 

0  y  p  0+2)  og 


(29) 


Z'  -  O+p)  (e’w  +  e  w')  -  fCAT'w  +  6Q'  (30) 

P  P  P  P  P  P  P 


The  non homogeneous  wave  equation  (27)  is  subject 
to  the  boundary  condition  of  the  form 


n  *  V  p ’  =  -f 


(31) 


where 


f  “  C  [  +  7^—  +  V(u*u*  +  u-u*) 

1 1  'it  -  -  -  ■ 


(u  x  '  *  +  u  x  '  +  u,x*')]*  n  -  H'*n-S’.n  (32) 


To  remove  the  time  dependent  terms,  we  introduce 

exponential  oscillations  of  the  form 

i  -  iakt 
p  =  p  e 


,  »  i ak  t 

u  =  u  e 


u*  -  \>*e 


i.ikt 


(3).i) 

(33b) 

(33c) 


(Q’  ,  ,  Z')  =  (Q,;,Z)e 


iakt 


(34) 


where  k  is  the  wave  number  defined  as 


-Z-  (cj  -  ia) 


(35) 


Substituting  (33)  into  (27)  and  (31),  respectively, 
we  obtain 


V2p  +  k2p  =  fi 
n-  Vp  = 


(36a) 

(36b) 


The  natural  modes  and  frequencies  corresponding 
to  (36)  are  determined  from 


V  PN  +  k  pN  =  ° 


n  •  VpN  =  0 


(37a) 

(37b) 


By  means  of  the  Green's  function  Integral 
and  using  (3  5)  -  (37),  it  can  be  shown  that 


=  kN  +  k  (  +  f  V0 


(38) 


where 


en  =  J  K' 


AZ 


(39) 


fi  =  i|  (0-7)^+  i  pN?  -u 


P  [  7 2 ( u • u '  +  u-u*)  -  7-  (ux?*+  u*xt 


+  u ’ x^ )  ]  +V.S’  +  7 •  H '  -  i^(Q'  +  i'+Z')(40) 


f  =  p [ iak(u '+u*)+7(u*u  +u- u) 
.  c  A 


(ux£*  +  Ox£  +  u'xf,)]-n  -  S*  »n  -H*n  (41) 


From  the  definition  given  by  (35)  and  (38)  and 
equating  the  imaginary  parts  it  follows  that  the 
growth  rate  a  takes  the  form 


a  =  nt  +  a 
P  g 


(42) 


where  rap  is  the  growth  rate  due  to  particle  damp¬ 
ing  and  refers  to  the  growth  rate  due  to  acous¬ 
tic  and  vortical  oscillations  of  the  gas  alone, 


and 


v-%1  t/((?'+?,)(I)-v) 

+  cj  (Q'+  $'+z')(R)pNJ;:  1 


(43) 


,  =  '  2K71  ^JK1'  +°*)(R>Pn  +  ?dF 
+  Y  J  P^G".,dr+*k“J  u  ((u'+u*)(I).  vJpjj.nd 

I  -  I 

-(2'i  +  l)J  PN(u'  ")PNJ  "''J  V  -“P^' 

“k^J  “■(<“>*>  (I->  )  PNd' 

‘  tv!  (("V/  +  +  U'xi)(1)-'.)p.. 


d  ]  (44) 


where'  (I)  end  (R)  indicate  imaginary  and  real 
parts,  respectively. 


Note  that,  in  arriving  at  the  stability  inte- 
trals,  it  is  necessary  to  perform  integration  by 
parts  "twice"  on  all  terms  associated  with 
the  convection,  viscous,  and  drag  forces  originated 
from  the  momentum  equation  as  the  wave  equation 
was  obtained  by  first  differentiating  the 
momentum  equation.  This  will  allow  appropriate 
Neumann  boundary  conditions  appearing  on  the 
boundary  surfaces.  A  glance  at  (43)  and  (44), 
however,  reveals  that  only  the  first  integration 
by  parts  was  performed  on  all  terms  other  than 
the  convective  terms  arising  from  the  momentum 
equations.  The  remaining  terms  are  yet  to  be 
integrated  by  parts  "one  more  time"  when  their 
forms  of  derivatives  are  explicitly  shown. 

3.  Evaluations  of  Particle  Damping  Stability  Inte¬ 
gral 


3.1  Approximate  Analytical  Solution 

We  now  return  to  (43)  and  examine  (16)  for 
perturbation  quantities.  The  viscous  drag 
forces  are 


where  T  and T  are  the  dynamic  and  thermal 
ation  times,  respectively, 

P  d! 

T  =  — - - 

d  IRu 


1  £_ 

2  C_ 


Prr 


relax- 


(55) 

(56) 


with  Ps  =  density  of  the  condensed  material  and 
and  d  =  average  diameter  of  the  particles. 


To  determine  the  effect  of  viscous  and  heat 
transfer  losses  arising  from  (48)  and  (49)  we  write 


-  1 


kN 


[  (S,(I)-?)pNd:>+r, 


JqV^i 


(57a) 


assuming  an  incompressible  one-dimensional  flow, 
we  arrive  at  an  analytical  solution. 


+r^  1  ( 5  71 

p  d  V  2p  Vpi- 

The  analytical  integraion  of  ("43)  is  not  pos¬ 
sible  and  it  is  our  approach  to  provide  numerical 
integration  via  finite  elements.  In  what  follows 
we  outline  the  detailed  procedure. 


H'  =  6f'  +  o' 

-p 


(45)  3.2  Finite  Element  Formulation 


where 

f.T'  =  -  p  (iak(u'  -  u’)  +  (u  •?)(;'  +  (u'-V.)uD 


-P  P  ~P 

-  (u  -7)0'  -  (G’-7)u] 


O’  =  (u  -  u)  w '  +  ( u '  -  u')w 

-  -p  -  p  ~p  -  p 

Simi  larly , 


(46) 

(47) 


S'  =  ii[V2u'  +  i  V(7  •  u')] 


(48) 


The  linear  viscous  losses  due  to  conduction 
heat  transfer  are  contributed  from  Q', 


Q'  =  (49) 

The  heat  transfer  interaction  energy  Z*  is  given  by 

Zf  «  Y’  +  <Q’  (50) 

P 

where 


Y*  =  (l+.'Xe'w  +  e  w')  -  8c(T  -  T')w  -  c'.u 
PPPP  P  P  -  ~ 


-  S  *  ur 


(51) 


The  finite  element  formulation  for  the  stabil¬ 
ity  integral  due  to  gas  oscillations  alone  given 
by  (44)  was  reported  in  [11).  Thus,  our  emphasis 
In  this  paper  is  concentrated  on  the  stability 
integral  representing  only  the  particle  damping  as 
shown  by  (43) . 


Substituting  (45)  -  (52)  Into  (43),  we  obtain 

(58) 


a  =  a„  +  a„  +  a„  +  aA+a, 

p  pH  pS  pQ  pt  pZ 


where 
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(n*  u)  (^y)Pjj-  (n-u')  (u-V)pN)dr 

^  fp(G(o;. 

S  L  pK'p  ~p 


'v)-  VpN  +  “P  (Vv)‘Vpn 


-  u-  (u'*y)  Vp  -  u'(u,V)'Vp„  -  (u  -  u)w'p 
-  -  N  -  N  -p  p 


p'-N 


(u'  -  u')w  p  )d0 
~p  ~  p  N' 


=  -PCf  |  iak (T •  -  T')  +  (G  • 

P  P  P 

V)T'  +  (u'  ■ V ) T 

P  P  P 

-  a20 

(O’  -  u  ' ) • 7 

)  PNd  '_ 

(I) 

(59) 

-  (G  •  V ) T '  -  (Ci '  ‘  ) T  I 

(52) 

P  - 

P 

Retaining  only  1  Fp  of  (45),  'Q'  in  (30)  we  re¬ 
write  (43)  in  the  form 


pS 


"  2e£  [^:J('.FJI>‘  )pNd  +  :\  IO}> 


from  which  a  simple  formula  m.iv  be  derived  for  a 
linearized  one-dimensional  problem  [5], 

-  I 


■-%  [4  J(<^ 'S 

r 

+  (n-V)V  •  (u*+u')pN)  d” 
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kN  W 


1) 


,+1"  1+('Vd* 


^  r  +  C7-1),P-  (54) 


'CP  1+<nV* 
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The  stability  Integral  n  Is  simplified  to 

P 


tin.  D(u*  +  u*) 

_ L _ 1 _ L_ 

3x  3x. 

j  i 


^  1  |  :iV"p  ‘  VPNd- 


'•u.  3(u|  +  u*)  (i) 

-23^— k - 


!pn  an 


=  7^  I  1(1+3)  (e'w  +  e  i')  -  BC(T'-T')w 
2tjl  J,  P  P  P  P  P  P 

it 
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where  the  superscript  (I)  denotes  the  imaginary 
part.  Note  that  the  boundary  terms  in  oipjjt  ipg, 
and  ipQ  have  appeared  in  addition  to  numerous 
other  terms  representing  the  viscous  forces  and 
thermal  dissipation  interacting  with  particles. 
Otherwise,  the  results  are  the  same  as  in  Culick[5], 

To  evaluate  numerically  these  integrals,  we 
introduce  three-dimensional  isoparametric  finite 
elements.  Any  variable  X  may  be  interpolated  as 


where  0^  is  the  interpolation  function  with  u  de¬ 
noting  the  global  node  number. 

The  procedure  for  the  solution  of  eigenvalues 
and  eigenvectors  for  both  acoustic  and  vortical 
motions  are  presented  in  (111.  Once  the  natural 
modes  and  frequencies  of  acoustic  and  vortical 
motions  are  determined,  then  the  final  form  of 
the  stability  integral  assumes  the  form 


i.: 


P  2E; 


+  I  cpcakNpp(Tp  -  r)pNdn  1 


In  terras  of  the  dynamic  and  thermal  relaxation 
times,  we  obtain 

-4 r-rV  »  — vf5» 
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d  pal^ 

1  *  ‘“T.  ,,  ten  s 


The  classical  acoustic  modes  are  of  the 
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a  *  m  mn 


k^  =  k2  +  k2 

N  i  ran 


where  m=0,  1,  2,  .  .  .  ,  k^  «  with 

J.  =  0,  1,  2,  .  .  .  and  k  are  the  roots  of  the 
derivative  of  Bessel  funcFion  J  (k  r)  such  that 


T-  [J  (k  r)  ] 
dr  m  mn 


0,  at  r 


r 

r  i  r  i 

.  -  -  f. 

G(;_)d')= 

*  o 

|  i  J  l  F(  -;,n)‘ir,dn 

P  ^  1  +  8  J 

.  .  2  _ 

Substituting  (68)  -  (72)  into  (67)  yields 


l+(uxd)  2  VPN 


G(  %*ti»  r  )  d  ~d  r  ,d  ry 


+  C(y-l)  Tt 


l+(.,»Tt>2 


+  —  p  lrdrdPdz 


Those  integrals  are  then  converted  to  Gaussian 
quadrature  of  the  form 

m  ro 

1  ’  Z  Z  MiWjf(ri.'G) 

i-l  i-1 

rn  m  m 

+  L  Z  £  wiwjV:(  i.g.V 

i *  1  j  =  l  k« 1 

win*  re  W.  represents  the  weights  and  f, .  ,  n-,  de¬ 
note  abscissa  of  the  Gaussian  quadrature  il2]. 

4 .  Example  Problems 

To  demonstrate  the  validity  of  the  procedure, 
we  consider  a  circular  cylinder  problem  (Fig.  1)  .. 

that  the  results  may  be  compared  with  the  analyti¬ 
cal  solution  given  hv  (*34). 


Numerical  analyses  have  been  carried  out  with 
the  following  assumptions:  £  *  0.2,  C(y-1)/C  =  1 

ps/p  =  1,  CPr/Cp=l,  Pr/(y-l)  =  1.  These  assBmp- 
tions  yield  Td  =  d2/18  and  it=d2/12.  The  results 
of  calculations  based  on  (54)  and  (73)  for  the  first 
axial  mode  with  various  frequencies  ( w  ”  300,800, 
1,800  Hz)  are  shown  in  Fig.  2.  Note  here  that  the 
trend  toward  decrease  in  stability  with  a  decrease 
in  frequencies  is  evident  for  both  one  and  two- 
dimensional  geometries.  However,  the  one-dimension¬ 
al  calculation  overestimates  (Fig. 3)  the  stability 
(particle  damping)  at  higher  frequencies  ( .-1 ,800  Hz). 
The  ranges  for  optimum  diameters  for  the  two-dimen¬ 
sional  cylindrical  system  (3  microns  for  gj  =1,800  Hz; 
4  microns  for  =  800  Hz;  6  microns  for  u  *  300  Hz) 
are  still  maintained  almost  the  same  as  in  the  one- 
d i mens i ona 1  approx i ma t i ons . 

Referring  to  Fig.  4  for  the  first  axial  and 
tangential  modes,  the  trends  are  s i gni f i cant  1 v  dif¬ 
ferent  f rom  these  of  the  axial  mode.  The  two-dimen¬ 
sional  cylindrical  system  is  more  stable  than  the 
onc-d imens ional  system  (Fig. 3)  for  all  frequencies. 
Although  the  ranges  of  optimum  diameters  remain  ap¬ 
proximately  the  same  as  in  the  rase  of  the  axial 
mode,  there  is  an  indication  that  they  shift  toward 


■  -j Mv! ‘.v: »" 


\  - -  <J  -  300  Hi 

\  _ =  800  Hi 


PARTICLE  DJA. 


(loglO  »  micron) 

Fig.  3  Growth  constants  vs  particle  diameters,  one-dimensional 
geometry 
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Oogjo  *  micron) 

Fig.  4  Growth  constants  vs  particle  diameters  for  first 

axial  and  tangential  mode  for  cylindrical  geometry 
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larger  sizes  as  the  frequencies  decrease. 

In  general.  It  Is  clear  that  one -dimensional 
apprnxmat ions  for  the  axisvmme trie  cylindrical  geo¬ 
metries  would  lead  to  erroneous  results.  The  trend 
of  deviation  depends  on  different  modes  of  oscilla- 
t ions. 

For  more  complete  analysis,  the  flow  field 
based  on  (1)  through  (22)  and  subsequent  evalua¬ 
tions  of  the  stability  integrals  (58-63)  must  be 
calculated.  Note  that  there  are  seven  variables: 
ui’  uPi>  T’  TPi’  pg*  pp,  and  p.  Therefore,  the 
required  governing  equations  consiste  of  (1),  (2), 
(9),  (12),  (17>,  and  (20).  The  drag  forces  are 
based  on  the  Stokes  law  given  by  Fp .  =  3*Tp  d(up^-Uj ) 
and  Qp*2  tt<  d(Tp-T)  where  d  denotes  the  average 
particle  d'ameter  and  <  is  the  thermal  conduct ivity . 
Initially,  all  incremental  quantities.  Sup.,  6Tp, 
^Fpi»  5Qp,  5ep,  and  6G  are  set  equal  to  zero. 

Then  these  values  are  calculated  from  the  initial 
solution  and  iterations  continue  until  convergence. 
The  results  of  this  analysis  will  be  reported  else¬ 
where  . 

5.  Concluding  Remarks 

A  rigorous  formulation  of  the  governing  equa¬ 
tions  involved  in  interactions  of  acoustic  oscil¬ 
lations  with  vortical  oscillations  have  been  demon¬ 
strated.  The  finite  element  analyses  for  both  flow 
field  calculations  and  stability  integral  evalua¬ 
tions  are  shown  to  be  convenient. 

Simple  example  problems  of  two-dimensional 
ax i symmetric  geometries  are  solved  and  compared 
with  one-dimensional  approximations.  It  is  shown 
that  there  is  a  trend  toward  decrease  in  stability 
with  a  decrease  in  frequencies.  However,  the  one¬ 
dimensional  calculation  overestimates  t lie  stability 
at  higher  frequencies.  The  optimum  ranges  of  diamet¬ 
ers  for  stab i 1 i ty ,  however ,  are  approximately  the 
same  in  both  cylindrical  and  one-dimens ional  geo¬ 
metries.  For  the  first  axial  and  tangential  modes, 
however,  the  trends  are  significantly  different  from 
those  of  the  axial  mode.  The  two-dire^s  i  ,inal  cylin¬ 
drical  system  is  more  stable  than  the  one-dimension¬ 
al  system  for  all  frequencies.  There  is  an  indica¬ 
tion  that  optimum  particle  diameters  shift  toward 
larger  sizes  a>  the  frequencies  decrease. 
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APPENDIX  3 


UNSTEADY  RESPONSE  OF  BURNING  SURFACE 
IN  SOLID  PROPELLANT  COMBUSTION 


T.J.  Chung*  and  P.K.  Kim** 
Department  of  Mechanical  Engineering 
The  University  of  Alabama  in  Huntsville 
Huntsville,  Alabama 


Abstract 

This  paper  is  an  update  to  the  earlier  paper 
concerning  the  multi-dimensional  calculations  of 
combustion  response  functions.  Our  ultimate  goal 
is  aimed  at  obtaining  more  precise  information  on 
combustion  responses  of  heterogeneous  solid  propel¬ 
lants  subjected  to  arbitrary  cavity  flow  fields, 
i.e.,  pressure-  and  velocity-coupling  with  turbu¬ 
lent  flames  and  possible  effects  of  radiation. 

This  paper  is  a  first  step  toward  such  a  goal.  If 
multi-dimensional  flow  fields  are  allowed,  the 
combustion  responses  can  no  longer  be  determined 
in  a  closed  form,  but  they  would  require  computer¬ 
ized  numerical  analyses  for  a  iarge  system  of  e- 
quations  as  a  result  of  first  and  second  order 
perturbations.  All  excited  frequencies  are  calcu¬ 
lated  by  means  of  eigenvalue  analyses,  and  the 
cumhustion  response  functions  corresponding  to 
these  frequencies  are  determined.  For  simplicity, 
the  calculations  are  limited  to  homogeneous  pro¬ 
pellants  and  laminar  flames.  Effects  of  radiation, 
incident  angles  of  impressed  pressure  waves,  and 
velocity-coupling  by  means  of  the  second  order 
perturbation  have  been  shown. 


Nomenclature 

a0  mean  local  speed  of  sound 

A  radiation  boundary  area  (Eq.  A-8) 

admittance  at  burning  surface 
B  frequency  factor 

c  specific  heat  at  constant  pressure  for  gas 

f.  activation  energy 

f  fuel  mass  fraction 

K  dimensional  frequency 

F  response  function 

h  heat  of  combustion  per  unit  mass  of  fuel 

H  radiation  function  (Eq.  A-8) 

I  radiation  intensity 

k  gas  thermal  conductivity 

K  dimensionless  wave  number 

•  flame  length 

I  latent  heat  of  vaporization,  radiation 

mean  free  path 
n  mass  flux 

Mach  number  at  burning  surface 
n  constant  exponent  ( Eq .  3,  A-6) ,  direction 

cosine 

N  conduction  to  radiation  parameter  (Eq.  5) 
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pressure 

Prandtl  number 

radiative  heat  flux  vector 

burning  rate 
gas  constant 

dimensionless  distance  (Eq.  16) 

Reynolds  number 

real  part  of  complex  variable  (Eqs.  33,34) 
time 

temperature 

gas  velocity  parallel  to  the  flame  surface 

gas  velocity  vector 

axial  core  velocity  (Eq.  16) 

gas  velocity  normal  to  the  flame  surface 

radiation  volumetric  element  (Eq.  A-8) 

reaction  rate 

coordinate  parallel  to  the  burning  surface 
coordinate  normal  to  burning  surface 
oxidizer-fuel  ratio 
thermal  diffusivity 
Eq.  (3) 

constant  exponent  (Eq.  A-6) 
perturbation  parameter 
Eq.  (7) 
eigenvalue 

solid  phase  eigen  function  (Eq.  B-17) 
Lagrange  multiplier 

dimensionless  radiation  source  function 

(Eq.  A-9) 

pressure  index 

Eq.  (A-14) 

mass  density 

optical  depth 

angle  between  normal  vector  to  the  surface 
(Eq.  A-8) 

finite  element  interpolation  function 

wave  incidence  angle 

eigen  vector 

frequency 

albedo 


Subscript  and  Superscript 


dimensional  quantity 
vector  quantity 
mean  or  constant  value 
solid  phase 
flame  edge 

radiative  surface  wall 
finite  element  global  node  number 
zeroth  order  perturbation 
first  order  perturbation 
second  order  perturbation 
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I.  Introduction 


2.  Governing  Equations 


The  response  function  as  related  to  the  com¬ 
bustion  of  solid  propellants  has  been  studied  by 
numerous  investigators.  Some  of  the  review  arti¬ 
cles  include  Culick  [1]  on  homogeneous  propellants 
[2-6]  and  Cohen  [8]  on  heterogeneous  propellants 
[9-14].  Both  of  them  are  limited  to  one-dimen¬ 
sional  and  quasi-steady  situations.  In  recent 
years,  some  works  toward  non-steady  problems  have 
been  attempted  [15-17],  The  majority  of  the  dis¬ 
cussions  of  response  functions  in  the  literature 
are  concerned  with  pressure-coupling  usually  ap¬ 
plicable  in  the  linear  stability.  Computations  of 
the  response  function  for  velocity-coupling  impor¬ 
tant  in  a  non-linear  process,  however,  remain  in  a 
state  of  infancy,  although  some  initial  attempts 
toward  this  subject  have  been  made  [19], 

As  noted  by  Flandro  [16],  the  velocity-coup¬ 
ling  may  be  accommodated  by  second-order  pertur¬ 
bations  of  mult i-d imensional  gas/solid  governing 
equations.  This  approach  deviates  drastically 
from  the  one-dimensional  analysis  which  has  been 
adopted  for  nearly  three  decades.  With  modern 
computers  and  various  tools  of  numerical  analy¬ 
sis  available,  however,  it  seems  possible  to  relax 
many  undesirable  restrictions.  Observations  indi¬ 
cate  that  combustion  oscillations  are  time-depen¬ 
dent  and  often  non-linear  as  influenced  by  turbu¬ 
lent  flow  fields,  which  may  lead  to  erosive  burn¬ 
ing  and  unstable  oscillations.  Unfortunately, 
however,  a  computational  tool  for  the  most  exact 
analysis  of  complicated  physical  phenomena  such  as 
these,  even  if  developed,  will  not  fit  into  the 
currently  available  computer.  Thus,  naturally, 
some  approximations  and  simplifications  must  still 
be  Introduced  to  any  theoretical  formulations  con¬ 
ducive  to  numerical  analysis. 

With  this  in  mind,  some  additional  materials 
and  reassessmen  s  concerning  the  multi-dimensional 
calculations  of  ombustion  response  functions  are 
presented  herein  as  an  extension  to  the  previous 
paper  [17].  The  influence  of  heterogeneity,  sur¬ 
face  roughness,  particulate  matter,  or  turbulence 
will  not  be  considered  at  this  time.  We  include 
the  effect  of  radiation,  adopt  a  simple,  premixed, 
single-step  laminar  flame,  expand  the  solid-gas 
governing  equations  into  first  and  second  order 
perturbations,  and  finally  perform  eigenvalue 
analyses.  Complicated  boundary  conditions  at  the 
solid-gas  interface  and  flame  edges  are  imposed  i- 
deally  by  means  of  Lagrange  multipliers  together 
with  finite  elements.  Calculations  are  carried 
out  for  various  incident  angles  of  impressed  pres¬ 
sure  waves  at  the  flame  edge  boundaries  [16]. 

Computed  results  show  that  the  natural  fre¬ 
quencies  are  clustered  around  low  frequency  range 
(^<10).  Spatial  distributions  of  field  variables 
corresponding  to  the  computed  frequencies  indicate 
that  the  oscillatory  behavior  is  pronounced  at  up¬ 
stream  and  gradually  diminishes  toward  downstream. 
Response  functions  oscillate  in  the  axial  direction 
with  peaks  occurring  at  the  midstream,  but  dimin¬ 
ishing  toward  flame  edges  for  both  first  and  second 
order  perturbations.  It  Is  also  shown  that  two-di¬ 
mensional  response  functions  are  multi-peaked  and 
may  become  negative  as  energy  sinks.  Finally,  It 
is  seen  that  the  effect  of  radiation  is  more  pro¬ 
nounced  In  the  second  order  perturbation  compared 
to  the  f  irst  order . 


The  combustion  of  solid  propellants  is  approx¬ 
imated  by  the  Arrhenius  law  and  a  one-step  forward 
chemical  reaction.  The  flame  of  erosive  burning 
is  assumed  to  be  a  premixed  laminar  deflagration 
with  a  single  species  (fuel)  [15,16]. 

The  basic  governing  equations  of  the  gas  phase 
are  composed  of  continuity,  momentum,  energy,  spe¬ 
cies,  and  state.  To  non-dimens ionalize  these  equa¬ 
tions,  we  proceed  as  follows.  Define  the  flame 
length  as 

e*  =■  k*/c*c*v*  (i) 

o  p  0 


where  *  denotes  dimensional  quantities,  subscript 
o  indicates  the  mean  value  in  the  chamber,  k  is  the 
thermal  conductivity,  p  is  the  density,  Cp  is  the 
specific  heat  at  constant  pressure,  and  v  is  the 
gas  speed  normal  to  the  surface. 

Introduce,  then,  the  following  dimensionless 
quantities: 


p  -  P*/P*  , 

p  - 

p*/p*  , 

T  •  t*/T* 

“  -  • 

t  • 

t*w*/l* 

o 

,  -  xj/l' 

“b  -  vX  . 

h  • 

h*/c*T* 

p  <x> 

,  E  -  E*/RT 

where  P  is  the  pressure,  T  Is  the  temperature,  u  is 
the  mean  flow  velocity,  t  is  the  time,  x  is  the 
length,  subscript  1  denotes  vector  quantity,  is 
the  Mach  number  at  the  burning  surface,  a£  is  the 

speed  of  sound  (a*  ”YfP*/ o*  with  y  •  c*/c*),  c  la 
the  specific  heat  at  constant  volume,  hpis  the  neat 
of  combustion  per  unit  mass  of  fuel,  subscript  00 
denotes  the  flame  edge,  E  is  the  activation  energy, 
and  R  is  the  gas  constant.  From  these,  the  govern¬ 
ing  equations  for  the  gas  phase  are  non-dimension¬ 
al  ized  and  the  explicit  forms  are  represented  as  in 
Appendix  A  (Eqs.  A-l  through  A-6) .  The  dimension¬ 
less  reaction  rate  explicitly  involved  in  the  equa¬ 
tions  of  energy  and  species  conservation  is  expres¬ 
sed  as  in  Eq.  (A-6).  Note  that  the  dimensionless 
frequency  factor  B  is  given  by 


B*k*T*"p*n 


B  - 


with  8 


* 

Ds 

3 


(3) 


where  8  Is  the  ratio  of  solid  to  gas  density,  n  is 
the  constant  exponent,  m  is  the  mass  flux,  W*  is 
the  molecular  weight  of  gas,  z  is  the  oxldlzer-fuel 
ratio,  *  is  the  constant  exponent  and  f  is  the  fuel 
mass  fraction.  It  is  noted  that  Eq.  (A-5)  is  valid 
under  the  assumption  that  the  perfect  gas  law  for 
the  reference  state  holds, 

P*  -  P  RT*  (4) 

o  o  o 


The  burning  of  solid  propellants  involves  ra¬ 
diative  heat  transfer  to  some  extent  as  the  combus¬ 
tion  chamber  is  an  emitting,  absorbing,  and  scat¬ 
tering  media.  It  may  be  assumed,  that  the  boundary 
surfaces  are  gray  and  diffuse.  Thus,  the  energy 
equation  is  given  by 


(u-T)T] 


Y-l 

Y  '<t 


+ 


i-n 

N 


(T* 


4tt 


H)  -  wh 


(5) 


where  u  is  the  albedo,  N  is  the  conduct lon-to-radi- 
ation  parameter,  H  is  the  radiation  function,  and 
n  is  the  dimensionless  radiation  source  function, 
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as  defined  in  Eqs.  (A-7)  through  (A-9) . 

For  the  solid  phase,  we  assume  that  the  mean 
burning  rate  is  constant  and  the  mean  surface  re¬ 
gression  is  given  by 

r*  -  P*V*/fl*  (&) 

O  o  3 


where  the  subscript  s  indicates  the  solid  phase. 
These  and  other  variables  are  related  as  follows: 


T 

s 


T  /T 
s  o 


*  * 

r  /v 


y 


* 

y 


*  *  *  /T 

p  c  V  /k 
o  p  o 


* 

o 


*  *  *  *  *2 

t  P  c  V  /k 
o  p  o  o 


*  *  ,  *  * 
k  c  /k  c 
s  p  OS 


(7) 


The  dimensionless  energy  equation  for  the  solid 
phase  takes  the  form  as  in  Eq.  (A-10) . 

Conversion  of  solid  to  gas  at  the  solid-gas 
interface  may  be  governed  by  an  Arrhenius  law  for 
the  dimensional  mass  flux 


* 

m 


B  P 


*nc 


(8) 


s  s 

with  E  ■  E*/RT*  being  the  dimensionless  surface 
s  s  s  ° 

activation  energy.  The  dimensionless  mass,  momen¬ 
tum,  and  energy  balances  across  the  interface  are 
given  in  Appendix  |A-11  through  A-13L 


3.  Perturbation  Expansion 

As  proposed  by  Flandro  [16],  the  perturbation 
expansion  of  all  variables  up  to  and  including  the 
second-order  would  enable  the  velocity-coupling,  as 
well  as  the  pressure-coupling,  to  be  adequately  mo¬ 
deled. 


Zeroth-Order  System 

All  perturbation  solutions  begin  with  the  ze- 
roth-order  system,  which  is  one-dimensional  and 
steady-state.  Thus,  we  have 
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where 
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For  simplicity,  we  may  investigate  an  effect 
of  an  impressed  pressure  wave  approaching  the  sur¬ 
face  at  arbitrary  incidence  given  by 


Solving  (17)  through  (22)  analytically  results 
in 
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and  the  velocity  distribution  is  given  by  (16). 
The  temperature  boundary  conditions  for  the  gas 
phase  are  obtained  as  follows: 

at  the  flame  edge  (y»°°) 


P  «  1  +  Ce^^fcosKfx*  +  x*)cos0  +  cosKysinO] 

+  € 2e^^Jt  [cosK(x*  +  x*)cos6  +  cosKysinQ ]  +  . . . 

(15) 

where  K  is  the  dimensionless  wave  number  (K  * 
2’TiF/a0vo),  u  is  the  dimensionless  frequency  (w  » 

2^'iF/v*2),  :  is  the  thermal  diffuslvity,  F  is  the 

dlmens?onless  frequency  in  Hertz,  0  is  the  arbi¬ 
trary  incidence  angle,  and  x*  is  the  antinode.  It 
is  also  proposed  (16)  that  a  simple  analytical  mo¬ 
del  for  the  mean  flow  streamline  in  the  vicinity  of 
the  surface  be  given  by 

— y/R 

u  *  uc(l  “  e  c)i  +  vj  (16) 

where  u  represents  the  axial  core  velocity,  and  Rc 
is  the  dimensionless  distance  from  the  surface. 
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at  the  solid-gas  Interface 
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Note  that  the  temperature  is  continuous  but  its 
gradient  is  discontinuous  at  the  interface. 


Perturbation  Systems 

Perturbation  expansions  are  performed  by  sub¬ 
stituting  Eqa.f9-14)  into  the  governing  equations 
using  the  following  perturbed  variables  and  Taylor 


series  expansion  about  the  origin  for  the  reaction 
rate  and  radiative  terms. 
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The  results  are  listed  in  Eqs .  (B-l  through  B-9) . 

The  boundary  conditions  corresponding  to  these 
perturbation  equations  include  the  flame  edge  (y*20), 
solid-gas  interface  (y*0),  and  the  location  deep  in 
the  solid  (y*-*).  At  the  flame  edge,  the  mass 
fraction  must  vanish  (Eq.  B-10) ,  and  the  impressed 
pressure  wave,  together  with  the  equation  of  state 
(Eq.  5),  will  be  specified.  The  temperature  bound¬ 
ary  condition  (Eq.  8-12)  arises  from  the  assumption 
that  the  flow  is  isentropic  at  the  flame  edge.  On 
the  other  hand,  the  mass  flux  fractions  of  the  fuel 
and  oxidizer  at  the  solid-gas  interface  are  assumed 
to  be  fixed  by  the  composition  of  the  propellants 
( Eq .  B-15) .  It  also  follows  from  Eq .  (A-13)  and 
Eq.  (B-9)  that  the  temperature  boundary  conditions 
at  the  solid-gas  interface  take  an  explicit  form 
as  shown  in  Eq.  ( B  —  1 6)  .  Furthermore,  as  a  conse¬ 
quence  of  mass  and  momentum  balances  across  the  in¬ 
terface,  we  arrive  at  the  velocity  boundary  condi¬ 
tions  (Eqs.  B-20,  B-21)  .  For  the  location  deep  in 
the  solid  (y»— *.)  ,  the  temperature  assumes  a  cons¬ 
tant,  independent  of  time,  as  defined  in  Eq.  (B-24). 

The  second  order  perturbation  system  and  the 
corresponding  boundary  conditions  can  be  derived 
in  a  similar  manner,  and  the  results  are  listed  in 
Eq .  C-l  through  C-9  and  F.q .  C-10  through  C-25, 
respec  t i ve 1 y . 

Finally,  the  response  functions  for  the  first 
and  second  order  systems  are  calculated  from 
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where  A^  ^  and  A^  ^  denote  the  admittance  of  the 
first  and  second  order  systems,  respectively. 


In  this  formulation,  the  governing  equations  (Eqs. 
B-l  through  B-8)  are  used  for  the  first  order  sys¬ 
tem  and  Eqs.  C-l  through  C-8  are  used  for  the  se¬ 
cond  order  system.  Note  that  the  pressure  degree 
of  freedom  is  excluded  from  computational  processes 
through  the  equation  of  state.  All  boundary  condi¬ 
tions  (Eqs.  B-10  through  B-22  for  the  first  order 
system;  Eqs.  C-10  through  C-25  for  the  second  order 
system)  are  cast  in  the  form  of  the  boundary  matrix 
equation , 
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where  r*l,2,...m,  m  being  the  total  number  of  bound¬ 
ary  conditions  (m=5  in  this  case).  Introducing  the 
concept  of  Lagrange  multipliers  Ar  and  minimization 
of  variational  principles,  we  now  arrive  at  the 
finite  element  matrix  equations  with  all  boundary 
conditions  properly  imposed. 


(38) 


This  represents  a  system  of  complex  characteristic 
equations  lending  itself  to  eigenvalue  problems. 
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These  expressions  (Eq.  38  and  Eq.  39),  of  course, 
represent  either  the  first  or  the  second  order  per¬ 
turbation  system.  Explicit  forms  of  these  equa¬ 
tions  are  shown  in  Ref.  (17).  Complex  eigenvalues 
and  eigenvectors  can  be  calculated  from  Eq.  (39). 
However,  the  system  given  by  Eq.  (38)  represents 
forced  oscillations  with  F,.  and  br  serving  as  forc¬ 
ing  functions.  Thus,  the  actual  amplitudes  corre¬ 
sponding  to  each  eigenvalue  (natural  frequency) 
can  then  be  calculated  from  Eq.  (38).  From  this 
information,  it  is  now  possible  to  determine  the 
response  functions  corresponding  to  the  first  and 
second  order  perturbation  systems.  Finite  element 
equations  are  summarized  in  Appendices  D,  E,  and  F. 


5.  Discussions 


4.  Finite  Element  Calculations 

Fig.  1  shows  the  domain  of  study  of  the  flame 
zone  (0  y  *)  with  the  boundary  conditions  to  be 
specified  at  the  flame  edge  (y-r)  and  at  the  solid- 
interface  (y«0).  The  solid  phase  (0'ys-«i)  is  e- 
llmlnated  from  the  calculation  by  providing  the 
sol  id - inter face  boundary  conditions. 

The  finite  element  equations,  by  means  of  i- 
soparametrlc  four-node  linear  elements,  are  con¬ 
structed  in  the  form 

( i-Ar^  ♦  Br, )X^  *  F -  (35) 

where  denotes  the  solution  vector  consisting  of 
nodal  values  of  the  density,  velocity,  temperature, 
mass  fraction,  and  radiation  function,  l.e.. 


Additional  calculations  beyond  the  previous 
paper  [17]  are  carried  out.  In  particular,  first 
order  eigenvalue  solutions  and  amplitudes  corre¬ 
sponding  to  each  of  the  natural  frequencies  are 
presented.  The  eigenvalue  calculations  for  the 
second  order  perturbation  system  are  not  included 
at  this  time,  although  the  amplitudes  for  the  se¬ 
cond  order  system  corresponding  to  arbitrarily  se¬ 
lected  frequencies  are  computed.  The  computation¬ 
al  domain  is  as  shown  in  Fig.  1,  and  the  various 
constants  used  are  as  follows:  E-10,  Eg-4 ,  L-0.15^ 

Ts-0.35,  T^0)-l.O,  Mb-0.01,  B-1000,  z-1,  C-Ol, 
"5*0.73,  k-1,  Y-1.27,  x*=0 ,  uc-l,  Rc-10,  Pr-1.0. 

The  results  of  the  eigenvalue  analysis,  as 
shown  in  Fig.  2,  represent  135  frequencies  (corre¬ 
sponding  to  all  variables  over  135  finite  element 
nodes).  Note  that  most  of  these  frequencies  are 
clustered  around  a  low  frequency  range  (oKlO.O) 
with  peaks  occurring  at  approximately  and 
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10' 3 .  This  trend  appears  to  be  true  for  both  rad¬ 
iative  and  non-radiative  cases,  although  cluster¬ 
ing  peaks  occur  at  slightly  different  frequencies. 

As  reported  earlier  [17],  the  spatial  distri¬ 
butions  of  variables  corresponding  to  the  comput¬ 
ed  frequencies  indicate  that  the  oscillatory  be¬ 
havior  is  pronounced  at  upstream  and  gradually  di¬ 
minishes  toward  downstream. 

Spatial  distributions  of  the  calculated  re¬ 
sponse  functions  for  the  first  and  second  order 
systems  are  shown,  respectively,  in  Figs.  3  and  4 
(ta-1.83,  >0.50,  N=*0. 1 ,  6-7T/2  and  Mb=0.01).  Al¬ 
though  it  is  premature  to  make  any  reasonable  as¬ 
sessments,  response  functions  tend  to  oscillate  in 
the  axial  direction  with  peaks  occurring  at  mid¬ 
stream  but  diminishing  toward  flame  edges.  This 
trend  appears  to  hold  true  for  both  first  and  se¬ 
cond  order  perturbations.  It  is  certain  that  this 
phenomenon  is  due  to  the  presence  of  mean  flow. 

Fig.  5  shows  the  frequency  dependence  of  re¬ 
sponse  function  for  the  first  and  second  order 
systems  at  location  M.  It  is  clear  that  two-di¬ 
mensional  response  functions  are  multi-peaked  and 
may  become  negative  as  observed  by  T'ien  (15]. 

The  notion  of  energy  sink  in  connection  with  nega¬ 
tive  response  functions  was  suggested  by  Flandro 
(lb],  and  this  may  perhaps  be  possible  for  an  os¬ 
cillatory  system.  The  effect  of  radiation  is  more 
pronounced  for  the  second  order  perturbation  com¬ 
pared  to  the  first  order.  In  general,  damping  ap¬ 
pears  to  prevail  at  low  frequencies  for  the  first 
order  system.  This  trend  is  reversed  for  high 
frequencies.  Large  negative  peaks  develop  for  the 
second  order  system  due  to  radiation  at  very  low 
and  very  high  frequencies,  with  intermediate  ran¬ 
ges  being  subdued. 

Effects  of  albedoes  on  response  functions  are 
shown  in  Fig.  6.  It  is  clear  that  response  func¬ 
tions  increase  with  a  decrease  of  albedoes  for  the 
first  order  response,  whereas  this  trend  is  not 
obvious  for  the  second  order  response.  On  the  o- 
ther  hand,  for  low  frequencies,  this  behavior  is 
reversed  for  the  case  of  the  first  order  response 
(Fig.  7).  Somewhat  erratic  behavior  is  observed 
for  the  second  order  response,  although  it  in¬ 
creases  negatively  with  a  decrease  of  albedoes. 

In  Fig.  8,  a  directional  dependency  of  res¬ 
ponse  functions  for  the  first  order  perturbation 
system  (at  location  M)  is  shown  in  polar  coordi¬ 
nates.  As  the  incidence  angle  increases  counter¬ 
clockwise  ,  response  functions  decrease  exponential¬ 
ly  to  the  minimum  at  '3«tt/2,  and  start  to  increase 
again  in  the  second  quadrant.  It  may  be  reasoned 
that  the  wave  incidence  parallel  to  the  burning 
surface  may  enhance  the  combustion  response  as  op¬ 
posed  to  the  case  of  normal  incidence.  Further¬ 
more,  Lue  influence  of  incidence  directed  upstream 
appears  to  be  stronger  than  the  incidences  direct¬ 
ed  downstream.  Similar  trends  prevail  for  the 
case  of  a  second  order  perturbation  system  (Fig. 9), 
except  that  the  response  diminishes  quickly  as  this 
angle  approaches  r‘* ^ . 

6.  Conclusions 

Eigenvalue  solutions  to  the  first  order  pertur¬ 
bation  system  associated  with  rocket  propellant 
combustion  have  been  carried  out.  Based  on  this 
information,  amplitudes  of  field  variables  in  the 
flame  zone  are  determined  and  response  functions 
for  selected  frequencies  are  computed  for  both 
first  and  second  order  perturbation  systems. 


The  following  conclusions  are  summarized: 

(a)  Most  of  the  natural  frequencies  are  clus¬ 
tered  around  low  frequency  range  (co<10.0). 

(b)  Spatial  distributions  of  field  variables 
corresponding  to  the  computed  frequencies 
indicate  that  the  oscillatory  behavior  is 
pronounced  at  upstream  and  gradually  di¬ 
minishes  toward  downstream. 

(c)  Response  functions  oscillate  in  the  axial 
direction  with  peaks  occurring  at  mid¬ 
stream,  but  diminishing  toward  flame  ed¬ 
ges  for  both  first  and  second  order  per¬ 
turbations  . 

(d)  Two-dimensional  response  functions  are 
multi-peaked  and  may  become  negative  as 
energy  sinks. 

(e)  The  effect  of  radiation  is  more  pronounc¬ 
ed  in  the  second  order  perturbation  com¬ 
pared  to  the  first  order. 

(f)  Response  functions  are  smallest  at  a  nor¬ 
mal  wave  incidence  and  increase  toward 
parallel  to  the  surface. 
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Appendix  A  -  Governing  Equations 


(Gas  Phase) 

Continuity 
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Species  Conservation 
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Appendix  B  -  First  Order  Perturbation  System 


(1)  Governing  Equations 
Continuity 

.  :(>)  .  .  "(0  (o). 


p<»  .  p(»)t(,)  +  a(i)T(0) 


Reaction  Rate 


~(i)  (o)  ~(  l )  j(l) 

•  [Jfsyr*"  -27tt  +277)i 


Radiative  Heat  Transfer 

v3(r!  “it  [«<t(0>>5*(i)  (B-7) 


H(l)  -  f  [4(l-n)(T(0))3  T(0 

h 

+  £»(,)l^dv 


Solid  Phase  Energ 


1  c  Ts  c  3y  5  3y 

_  7jt^  »  0 


(2)  Boundary  Conditions 


(A-ll) 

Flame  Edge  (y-00) 

(A-12) 

f(1)(oo)  .  o 

rL 

(A-13) 

p(0)(»)T(l)(«>)  +  0 

■  cosK(x  +  : 

l  -  k*/k* 

(A-14) 

o 

#<’>  -ii£_  „* 

v  '  y  Ui  oy  j  y-  It* 

u^'\")* - —  sinK(x*  +  x*)cos6 

fMb  0 


sinK£*sin9 


(B-10) 


(B-12) 


Solid-Gas  Interface  (y»0) 

,;(D  •>.<»>  -(i 

f(1)  .  (  ii -  )  +  (  M _  )  r  v_ 

1  +  '  3y  ;+  v  3y  V  T 


V  +  -^bx,-  L] 


v  5y  jf  "b‘  gcj"  go1"'1  ‘ 

X,  .  E  EL  ,  E 

-f  (1+^-^-B)  +  —  -  -i — 5-A 

^  8oT  tb  6o7 


where  A,  -  [  1  +  (1  +  4180 1/2 ] 

T  -  T 

s  a> 


T  -  T 
s  « 


Energy 

+  1 

12op(‘>T(2)  +  p(0)[(u(0).V)T(2)+  (u(0 

+  (u(2) •V)T(°)]  +  lwp(l)T(l) 

(B-15) 

+  P(l)[(u(0)*V)T(l)  +  (u(1)-V)T(l,)] 

_ 

+  P(2>(u(°>  *7)7  +  12 u  P(2)  - 

+  v3u)  -w(2)h-o 

Species  Conservation 

(B-16) 

12u)p(o)f(2)  +  p(0)  ((u(0) -V)f(2)  +  (G(l)  ■ 

(B-17) 

+  (u(2)  *V)f(l>)]  +  lop^f*'5 

(B-18) 

+  P(1)((u(0)-V)f(l)  +  (u(l)*V)f(0)] 

(B-19) 

+  p(2)(u(0)-V)f(0)  -  V2f(2)  +  w(2) 

U(0+  ’  '  ~k  (=^f(l  +  +  n3XiV_VtCOSK(X' 


p(2)  .  pC){(*)  +  3(0*0)  +  sc*)T(o) 


+  x*)cos0  +  sin0) 

V(l).  -  { [  (E  +  l)T(l>  +  (n  -  DTI 
+  s  +  s  s 


(B-20) 


1  ,  Es  c(D 

-  i^(r  +  n 

A8 

+  x*)cos8  +  sln0] 


n  X^V=)McosK(* 


O 

(B-21) 


+T^1^  +  ,<  V'’  -  cosK(x*  +  x*)cos8  +  8ln8 

(B-22) 


Deep  In  Solid  (y--00) 
T(l)(-»)  .  0 


(B-23) 

(B-24) 


Appendix  C  -  Second  Order  Perturbation  System 

(1)  Governing  Equations 
Continuity 

12o3<2>  +  V.(p(°>u(2)  +  S(l)u0) 

+  o(2)u(0))  -  0  (C-l) 

Momen  turn 

12UC(5>u<2)  +  P(0)[(u(°).V)u(2)  +  (lI(l).V)u"<1) 

+  (u<2).V)u(0)i+  iJ'W"' 

+  C(l)[(u<0)-V)u(l)  +  (u(1)-V)u(0)) 

.  CO),  (0)  (0)  1  -(2) 


Reaction  Rate 

;(*)  .„(»>{  iij!L  +  nil.  [2p(,>.  2 

(TC°))2  (T(°>)2 

+  2  +  (P(0)2  +  2P(2)  +  (  )2 

£<*)  ;(0  ,  i(D 

*27(5T 

-O)  T0)  f(0  ET(I) 

*22<  ’'-'jw*2??*^1 

,  T(l)  rf(D  ,  f(l)  .  ET(1)  . 

- 2  70  [2P  +  2  70  +  ^7577. ] 


.  ,  f(1)  r,f(0  .  T(0  ,  ET(1)  . 

+  27Td-[2p  '  2  T(°)  +  ^To)“ 


*77^°’ 


*(>)  2(0 

2  +  2  -Vv  ] 

T(»)  f(»)  ' 


Radiative  Heat  Transfer 


V3(R)  "  '^T  I4(T<0))3T(2)  +  6(T(°V  (T(0)2 

(C— 7) 

H(2>  -  f  [4(1-D(T(0))3  T<2)  +  6(T(0))2(T(,))2 

'v 

+  ^”(2>1^'dV  (C-B) 


Solid  Phase  Energ 


(2)  Boundary  Condition 
Flame  Edge  (y**00) 

f(2)(»)  -  0  (C-10) 

d<°>(oo)T(2)(<»)  +  3(‘>(<»)t(1)(<»)  +  S(2)(»)t(o)(») 

-  cosK(x*  +  x*)cos0  +  cosK£*sin8  (C-ll) 
o 


_  i=i  p<*>  _  A  I 

y  “  3y  ly.n* 

u<2\«)  «  -  sinK(x*  +  x*)cos0 


(C-12) 


(C-13) 


v(2)(»)  -  -  sinKS*sin0  (C-14! 

™b 

Solid-Gas  Interface  (y*0) 

0(2)  a;(0  *(»). 

f“' )  [  vh!  (4^)2 


!  .  C  ==•  t  +  +  ns)(  >++  [2ns-r 


+  n  (n  -i)  +  (JLi<'>+)2  +Im<2Vk-^J 

8  s  t.  c 


p(o)?(2>  +  -(OjO)  +  s(2V°j 


cosK(x*  +  x*)  cos  0  +  sin0 


Deep  in  the  Solid 


T(2)(-°°)  -  0 

S 


3x  |y— oo  3y  |  y- 


Appendix  D  - 


Finite  Element  Equations  of  First  Order  System 


2(2)  2(>).  , 

T_++  (  1=-^)2 


ai(2>  1 

(  —T~ -  ),  +  -n-CG  ~ 


6u  (Buy 


-  D(G+n  )  -  LG 


(C-17) 


(C-18) 


(C-19) 


E  ,  E 

—  +  L—  -  -J-zp  C 
C  Tg  SUTs 

T  -  - 

where  C  »  -3  ^ -  (Xi  -  —  ) 


G  -  (  =2-T(l)+)2+  2n  — T(l) 

Ts  s  Ta  + 

+  n  (n  -  l) 

s  s 

+  x  ^  s  3y  + 


+  [2n  —  T(l),  +  n  (n  -  1)  +  (=#T(l)+)2 


E  a  (») 

+  =~  T( 2  +  )  (  ^ -  )  }[cosK(x*  +  x*)  cos 0 

x  Jy  +  ° 

+  sin9)  2  (C-20) 

v^2  +  «  (R  -  I)  [cosK(x*  +  x*)cosP  +  sin0]2 


(C-21) 


F.  (E  +1)  ,,, 

R  -  (E  +1)T(  5.  +  -~5 -  (T  +)2 

S  T  T 


+  [E  (2n  -  =k  )  +  (n  -  l) ]T(  1  \ 


(1)  Governing  Equations 
Continuity 

<**<*  +  +^Tub1)  ‘  0  (D 


(iuAaB  +  DaB)uBi  +  EaBi  PB  FaSi  TB 

(D 

Energy 

<lh*^  +  +  -  (iu5^  +  *WPB 

'  "b0  ■  °  (D 


Species 

(iuAaB  +  Qa8)f8  +  RaS  TB  +  SaB  PB 


Radiation  Function 

tb°  +  <RP«s  -  V’b0  ■  0  (D 

Note  that  explicit  forms  of  these  matrices  are 
shown  in  Appendix  F. 

(2)  Boundary  Conditions 
Flame  Edge 

f£!)  -  0  (D 

(  z,  ab  +  Vtb°  -  J  (D 

T<0)bbpB,)  +  p(0)BBTs,)  “  ®  (D 

BaUo  ;  -  iH  (D 


( 1)  - 
Vs  '  11 


where 


“3  3y 

J  1  V“ 


y-i* 


B3  *  Vy-S* 

A  ’  _  t^SBS 


B  * 

G1  ■  cosK(x*  +  x*)cos0  +  cosK£*sin0 
o 

H  -  -  ^  [sinK(x*  +  x*)cos8]gBg 

T  -  -  slnKJl*sine 
™b 


(0-10) 

(D-ll) 

(D-12) 

(0-13) 

(D-14) 

(0-15) 

(D-16) 

(0-17) 


Solid-Gas  Interface 

<eb-  V'e0  -  (Gr  >vb0  +  (?r  >tb° 

3  3  (D-18) 


{^#(Tb-a)eb+  t(T  +  r  LVV,T|(,) 


gfs  ■  4  --B  "  f.  Ts  B'  ■  B 

i  ns 

■  U  T  <A-®V  -  nsL  (0-19) 


E6uS 


u)  Sr 


6  6 


,  o  u 
i  3  C  ~ 


U) 


Momentum 

(%  +  y“Si  +  (Vi  + vu)pS.2) 
+  vT  TJJ)  ■  -iwFBa  -  FC„ 


(E-2) 


Energy 

2“  a 
Y  AaS 


(i  %  A„n  +  GC^g  +  HHag  +  RAag)  Tg2)  -  00^f<2) 


aB 


6 


-  RB^H^  -  (12“Sag  +  V)pJ2) 


-iuFD  +  FE 

a  a 


Species 

(120^  +  GG^g  +  00C(g)fJ2)  +  HHagTJ2> 
+  «aePS2>  ■  iuFGa  +  ™a 

Radiation  Function 

ifaBTB2>  +  ®aB  ‘  AaB>HB°  *  'FZa 


(E-3) 


(E-4) 


(E-5) 


Note  that  explicit  forms  of  these  matrices  are 
shown  In  Appendix  F. 

(2)  Boundary  Conditions 
Flame  Edge 


8 


(E-6) 


EBvb‘)  -  B[(Esn) 


.  E  T  -T 

1  _s  ,  _s _ « 

Boj  -j*  C 

A  8 


>]EaT, 


.  n  T  -T 

--H-V*d  +  (v1)is 

(T<,iv4°  +  (p(°lE8)T8l)  -  ° 


where 


8  ’  *3  ‘  y 

3t>0 

F  -  — -S I 
f8  3y 

3f 


y-o 

(o) 


( 


3y  )sEe 


D  -  ( G2 ]  Eg 


G2  «  cosK(x*  +  x*)cos8  +  slnA 
o 


(D-20) 

.O) 

(i  ^  ab  +  bb)tb2>  ’  A 

(E-7) 

s 

+  -  B 

(E-8) 

(D-21) 

bbub2>  ’  1B 

(E-9) 

( D— 22 ) 

bsvb2>  ’  i7 

Where 

Afl  - 

6  3y  ly-e* 

(E-10) 

(D-23) 

(E-ll) 

(D-24) 

BB  -  Vy-£* 

(E-12) 

(D-25) 

Brg-p(0)c»)  bb 

(E-13) 

(D-26) 

B2a  -  T(<>)(”)  Bg 

(E-14) 

(D-27) 

A  -  ^  [Gl]  .Bg 

(E-15) 

Appendix  E  - 

Finite  Element  Equations  of  Second  Order  System 

(1)  Governing  Equations 
Cont lnulty 


(12wA  ,  +  B  a) 


»S  8 


+  C 


—  (?) 
»£1  uBl 


FA 


(F.-1) 


B  •  f  G1  -  *.»  i  j  3o^ 


(E-16) 


Cl  -  cobK(x*  +  x*)c os0  +  cosK5*sinP  (E-17) 
o 

H  -  -  [ sinK( x*  +  x*)cosP)  ,Ba  (E-18) 

>  rt  o  i*  o 

sinK  >  *slnf*  (E-19) 


41 


Solid-Gas  Interface 

<Ee-Vf62)  -  G 

ESUS  )  "  E8T82>  "  “ 

esvb2)  +  Vb2)  *  lE  +  5 


(ps + 1  —  >n  -  is“ +  t 

hA2)  +  *Tbtb0  * « 


(E-21) 


(E-23) 


bl  -  (a  +  16BJ<jV)1/4sln  -f  )  (E-41) 

C,  -  cosK(x*  +  x*) cos 0  +  slnO  (E-42) 

j  o 

E  E 

c,  -  2n  J-  T(,)  +  n  (n  -l)  +  (^T(l))2  (e-43) 
3  8  t.  •  •  Ts 


Gj  ”  ~  (T(l))2  +  (-  =2j+  <n  -1)/T 


Note  that  0-  denotes  the  principal  angle  of  the  so¬ 
lid  phase  eigenfunction  X  . 


8  '  VB 


eib"  t(0)(0)eb 

0(o)(o)Eb 


(E-27) 


Matrices  for  Appendices  E  and  F 


BaB  +  CaS 


8  3y 


c6--r=i<cj,‘'i! 


4  T  -T  E 
n-  _  /  1  s  o"  i 


(E-28) 


°B  “  {  -  V"  f  '  <Es+1)  }(C3)2  E8  (E-30) 


P6.  {i^  +  ^  +  }Eb.Kb 


E„  T  -T 


°8  "  *'i(ai  _  q  5  +  bi}  gj  zZ  ES  (E"32) 


“  ,  3f0)  v(0  .  .  T<*>  ,  _ 

1  ~W~  "V  "fr  sEe  ( 

5-i<c3>2  R£'(rT<l)  +  vV 

c  4s  7 

+  C3]BE2  < 

T  -T 

+  c3  —  1  .E8  ( 

0  *  (C3)2[TsG3  +  W?  ( 

l  1  1 

s  ”  "SZi  (al  *  c  )  ^C3*3E°_  bl^r’3 

+  VaV  < 


(E-33) 


(E-34) 


(E-36) 


d(5)bo8+  p(°>CaB+J  PrDae 


Ea  81  +  yMb  CaBl 


YM^  SSI 

T<°> 

YH^  CaBl 

^A 

Y  aB 

°(°)ca8  +  Dcx8-h«^ 


111  T<°>A 

Y  AaB 

2hv(0)T(°V 


p(°>Ca8+  DaB 


~  T(») 


+  20  2  ]A 


2„(0>T(5)A 


(  °) 

P  CaS  +  D*B 


T  -T 

~7T~  'bl[C313E3  BSc  Ial  ‘  \  )  tC3 


+  r  1  E  !  .  t[C,]aEp 


H  -  (C, 


:( l)T( 1} )gEg 


*  -ft  fl  +  (1  +  16B2u.2C2)1/4cos  }  (E-40) 


4  1=2  t(0)SA 
N  Ani8 

ikd  A 
4  it  N  .B 


-  4(i-i;)T 


’z 


RDae  “  4W  2aS 

T-  .  2w(0)T(0)A 
sa6  w  aB 

AaS  “  *a*B  dfl 

BaS  ■  L 

■'il. 

C  Q  -  [  u$0)4>  4>  ,  df2 

ctB  i  a  8,i 

C  ,  -  f  4>  $a  .  dQ 

aBi  Jn  a  8,1 

D  a  -ft  ,  d« 

aB  Jn  a,j  6,j 

Ect6l  "  jn  ui,juj  J  Vb  dfl 

c  .Ill  t(0)a 
b™e  y  1  Aae 

T  ■  2hw^°^T^°^A 

aB  a8 

Za8  *  2l.  VY[RTTlYdn'L  *a  dfl 


fO>  2  T(l)  2 

*<7ra> 

To>  f(»)  T(>) 

.2I_[2P()  +  21_  +  E_J 

rn  T(l)  f(1)  T(1) 

+  2P  t_2  70  +  2  7W+  e7To^ 

,  f(l)  f  ,  T(l)  (1)  T(l)  , 

+  2  joy [-2  +  2?  +  et^ 

ET(0  ,  ,  T0)  (I)  ,  f(0  , 

+  ['2  77)  +  2P  +  2  1 

v.p(ov(1)t(:)  +  p(iv(,)t(;)u^w^2 

1  1  ,1  N 


z  -  p(0V(,)f(!>  +  pOM’W 


i  .i 


i  .i 


tt/2  -t/cOS<>2 

RTT‘  jo  ^ d*g 

2w(0) 

°°a6  ”  f(o)  AaB 

PA«"-t,IPH>U*,>+P<‘)u^]A*a  dfJ 

"a  -  l  ^(,>“i,>1YVadn 

-  4v«*° 'M3-S0 

+  ^0<1>T^Va  dn 

Ft)  -  -  f  [0(1)T(1)1  $  $  d« 

a  L  Y  Y  a 


Ea  -  |  [hw0)X  -  Y], 


$  $  dft 


'  $  d« 

Y  a 


H  -  -f  [w(0)X  +  Z]  <t  <t  d« 

a  -  ‘Y  Y  a 


FZr(  *  12T(0)2|  _  ^[RTTlyW'J  't'n[T(1)2]rl4>a  dfl 


X  -  E-T'.-  ■■  (2P(l)  -  2  +  2  — .-r-  ]  +2o(lJr(l) 

-  <T(0V  T<°>  f(s) 
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